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1.0  INTRODUCTION 


Ergodicity  is  the  condition  which  enables  time-averaged  statistics  of 
random  processes  to  approximate  those  obtained  by  ensemble  averages.  Although 
this  condition  is  often  assumed  in  estimation  and  other  signal  processing 
applications,  the  dependence  of  the  ergodic  behavior  of  random  processes  upon 
fundamental  process  characteristics  (such  as  the  temporal  and  cross-channel 
correlation,  and  the  process  variance)  and  its  implications  on  estimation 
performance  evaluation  is  not  often  considered.  The  ergodicity  condition  for 
auto-  and  cross-channel  correlation  functions  expressed  in  terms  of  these 
fundamental  process  characteristics  is  derived  in  this  paper.  Specifically,  analytic 
expressions  are  developed  for  the  variance  of  the  biased,  time-averaged 
correlation  functions  for  stationary  discrete  complex  processes.  If  these 
variances  approach  zero  in  the  limit  of  infinitely  large  sample  sizes,  the  ergodic 
condition  holds.  The  expressions  derived  here  pertain  to  the  general  case  of 
Gaussian  processes  with  unconstrained  quadrature  components  where  the 
bandpass  processes  are,  in  general,  non-stationary  [51.  In  addition,  the  analytic 
expressions  are  simplified  for  the  special  case  of  complex  processes  with 
constrained  correlation  between  the  Gaussian  quadrature  components  yielding 
stationary  narrowband  bandpass  processes. 

The  expressions  developed  here  for  the  time-averaged  correlation  function 
estimators  provide  a  performance  measure  which  can  be  used  to  specify  the 
window  size  of  the  observation  interval  required  to  achieve  a  specific  value  of 
this  variance.  As  an  example,  the  increase  in  sample  window  size  required  to 
achieve  a  specific  value  of  the  variance  of  the  time-averaged  correlation  functions 
is  quantitatively  related  to  the  increase  in  temporal  correlation  as  well  as  the 
variance  of  the  underlying  process.  Validity  of  the  analytic  expressions  is 
presented  using  a  multichannel  process  synthesis  method  described  in  [1,21. 

Next,  Monte-Carlo  simulations  are  used  to  measure  the  error  variance  of 
several  parameter  estimators  using  zero-mean  time  series  whose  covariances  are 
functions  of  a  parameter  vector.  Both  Gaussian  and  non-Gaussian  processes  are 
considered.  In  this  analysis,  the  parameter  vectors  are  varied  to  simulate 
processes  with  various  temporal  and  cross-channel  correlation;  ie.,  processes 
ranging  from  narrowband  to  broadband.  The  dependence  of  the  parameter 
estimators  upon  process  correlation  is  contrasted  with  that  of  the  time  averaged 
correlation  function  estimators.  It  is  noted  that  the  performance  of  these  two 
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classes  of  estimators  have  markedly  different  dependence  on  process  temporal 
correlation.  Finally,  the  performance  of  the  parameter  estimators  is  compared  to 
the  exact  Cramer-Rao  bound. 
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2.0  PROCESS  DEFINITION 


In  this  chapter,  we  define  the  multichannel  autoregressive  (AR)  process, 
the  relation  between  the  AR  coefficients  and  the  correlation  matrix  of  the 
processes  known  as  the  Yule-Walker  equation,  a  wide  class  of  non-Gaussian 
processes  known  as  Spherically  Invariant  Random  Processes  (SIRP’s);  and  finally, 
a  brief  description  of  the  procedure  used  to  synthesize  these  processes  enabling 
the  control  of  the  temporal  and  cross-channel  correlation. 

2.1  MULTICHANNEL  AUTOREGRESSIVE  (AR)  PROCESSES 

The  multichannel  forward  autoregressive  (AR)  process  x(n)  is  expressed  as 
M 

x(n)  =  -  X  A^(k)x(n-k)  +  u(n)  (2.1) 

k=l 


where  A”(k)  is  the  kth  JxJ  matrix  coefficient  for  the  Jxl  vector  AR  process  x(n) 
of  model  order  M  and  jyL(n)  is  the  Jxl  forward  white  noise  driving  vector  with 
covariance  matrix  [Sflu-  We  note  that  A^(k)  is  expressed  in  terms  of  the 

Hermitian  operation  for  notational  convenience,  but  is  not  restricted  to  be  a 
Hermitian  matrix.  It  will  be  shown  later  that  this  notation  enables  us  to  express 
the  Yule- Walker  equation  in  terms  of  row  vectors  containing  the  Hermitian 
transpose  operation.  For  the  single  channel  case,  the  terms  in  eq(2.1)  are  scalers. 
Thus, 


M 

x(n)  =  -  X  a*(k)x(n-k)  +  u(n)  (2.2) 

k=l 


The  vector  ii(n)  from  eq(2.1)  can  be  expressed  as 

ll(n)  =  CvV(n)  (2.3) 

where  v(n)  is  a  Jxl  white  noise  vector  with  JxJ  diagonal  covariance  matrix  D^. 
Eq(2.3)  enables  us  to  obtain 
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lEflj,  s  E[u(n)uH(n)l  =  ElC„v(n)vH(n)Cvl  (2.4a) 

=  C,,E(v(n)vH(n)|Cv  (2.4b) 

=  C^,D^,Cv  (2.4c) 

For  the  special  case 

Dy  =  E[v(n)vH(n)]  =  I  (2.5) 

where  I  is  the  JxJ  identity  matrix,  eq(2.4c)  reduces  to 

[If]„  =  C^,Cv.  (2.6) 


Eq.(2.6)  indicates  that  the  JxJ  matrix  Cy  can  be  obtained  by  the  Cholesky 
decomposition  of  The  above  relationships  are  utilized  in  the  process 

synthesis  procedure  described  in  section  2.4  and  reference  [2].  Alternatively, 
eq(2.3)  could  have  been  expressed  as 

u(n)  =  Lzz(n)  (2.7) 

where  is  a  unit  diagonal  lower  triangular  matrix  and  the  vector  process  z(n)  is 

a  white  noise  vector  uncorrelated  in  time,  whose  components  are  also 
uncorrelated  across  channels.  The  covariance  matrix  of  u(n)  is  expressed  as 

=  E[u(n)uH(n)]  =  E[L2;Z(n)zH(n)L2]  (2.8a) 

=  LzE[z(n)zH(n)]Lz  (2.8b) 

=  (2.8c) 

where 

Dz  =  E[z(n)zH(n)].  (2.9) 

Eq(2.8c)  denotes  the  LDL^  decomposition  of  f2^f]u  provided  is  a  real, 

diagonal  matrix.  Furthermore,  the  channel  variances  are  contained  along  the 
diagonal  elements  of  D^. 
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2.2.  THE  YULE- WALKER  EQUATIONS 


In  this  subsection,  we  present  the  relationship  between  the  matrix 
coefficients  A^(k)  of  eq(2.1),  the  covariance  matrix  of  the  forward  AR 

white  noise  driving  vector  defined  in  eqs(2.4),  and  the  known  correlation  matrix 
[R^jj]  of  the  vector  x(n).  The  latter  is  defined  as 

[RJ  =  E[x^^  U  ^  )H].  (2.10) 

where 

=  [x^(n-M)  x^(n-M+l)...x^(n)]  (2.11) 

~n-M,n 

We  first  introduce  the  reversed  order  correlation  matrix  [RxxIm+i  defined 
as 

=  E[i„  (2-12) 

where 

=  [x^(n)  x^(n-l)...x^(n-M)]  (2.13) 

n,n-M 

and  ~  denotes  time  order  reversal.  We  note  that  the  matrix  [Rxxl  is  the 
correlation  matrix  of  x  ..  =  x  in  contrast  to  [Rxx]  which  is  the  correlation 

— n-M,n  — n,n-M  — 

matrix  of  the  vector  „  defined  in  eq(2.1 1). 

The  relationship  between  A”(k),  [Sflu  and  [R^]  is  expressed  by  the 
augmented  Yule-Walker  equations  [9,10];  ie.. 


A"  [Rj^l  =  {[Zf]|;  (01.401)  (2.14) 

where 

a”  =  [I  A”(1)  A”(2)...A“(M)]  (2.15a) 

[^]u  =  E[il(n)u”(n)]  =  [%]„  (2.15b) 

ll'r(n)  =  [Ui(n)  U2(n)  ...Uj(n)l  (2.15c) 
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The  corresponding  equation  for  the  wide-sense  stationary,  backward  AR 
process  is  expressed  as 

M 

x(n)  =  -  X  BH(k)x(n+k)  +  Ub(n)  (2. 1 6) 

k=l 


where  we  specifically  denote  the  backward  white  noise  driving  vector  U|j(n)  with 
the  subscript  b.  The  covariance  matrix  of  ujj(n)  is  expressed  as 

[Zbl„  =  E[U|,(n)  Ub(n  - 1)  (2.17) 

The  corresponding  equation  for  the  backward  Yule-Walker  equations  are 
expressed  as 

B“[Rxx1  =  ([01...(01  [ZblS) 

where 

2.3  NON-GAUSSIAN  AUTOREGRESSIVE  PROCESSES 

In  this  section,  we  discuss  the  generalization  of  the  single  channel 
autoregressive  processes  defined  in  eq(2.2)  to  a  class  of  non-Gaussian  random 
processes  known  as  spherically  invariant  random  processes  (SIRP)  [7,16,17,18]. 
The  non-Gaussian  form  of  the  random  process  x(n)  is  introduced  through  the 
white  noise  driving  term  u(n).  Following  Rangaswami  [16],  we  first  define  a 
spherically  invariant  random  vector  (SIRV)  as  a  random  vector  (real  or  complex) 
whose  PDF  is  uniquely  determined  by  the  specification  of  a  mean  vector,  a 
covariance  matrix  and  a  characteristic  first  order  PDF.  A  spherically  invariant 
random  process  (SIRP)  is  a  random  process  (real  or  complex)  such  that  every 
random  vector  obtained  by  sampling  this  process  is  an  SIRV.  An  important 
theorem  in  the  theory  of  such  processes  is  the  representation  theorem  [Yao]  stated 
as  follows. 


(2.18) 

(2.19) 
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Theorem  1  If  a  random  vector  is  an  SIRV,  then  there  exists  a  non-negative 
random  variable  S  such  that  the  PDF  of  the  random  vector  conditioned  on  S  is  a 
multivariate  Gaussian  PDF. 

For  the  simulation  of  an  SIRV,  we  consider  the  product 

(2.20) 

where  Uj  =  [uj  U2  ...  u^]^  denotes  the  SIRV,  Zj  =  [z,  Z2  ...  Zj^J"^  is  a  Gaussian 

random  vector  with  zero  mean  and  covariance  matrix  M  and  S  is  a  real,  non¬ 
negative  random  variable  with  characteristic  PDF  f^Cs).  Statistical  independence 
between  Zj  and  S  is  assumed  for  convenience.  In  [7],  several  characterisic 
PDF’s  for  f5(s)  are  considered  which  provide  various  PDF's  for  fjj(u).  Among 
others,  they  include  the  Chi,  Weibull,  Generalized  Rayleigh,  Rician,  the  K- 
distribution,  Laplace,  Cauchy,  Student-t  and,  as  a  special  case,  the  Gaussian.  In 
section  5.3.2,  we  consider  K-distributed  processes  using  a  form  of  the  Gamma 
distribution  for  fs(s). 

2.4  PROCESS  SYNTHESIS  PROCEDURE 

In  this  section,  we  briefly  describe  the  method  used  to  synthesize  the 
random  processes  used  in  this  study.  The  procedure  utilizes  the  multichannel 
Yule-Walker  equation  expessed  by  eq.(2.14).  Essentially,  the  desired  temporal 

and  cross-channel  correlation  are  specified  in  the  covariance  matrix  [R^], 

14  H 

eq(2.14)  is  solved  for  the  vector  of  matrix  coefficients  A  and  the  matrix  [Zf]u. 
The  A^(k)  coefficients  contained  in  are  used  in  eq(2.1)  directly  for  the 
desired  order  M  of  the  process.  The  matrix  [Zfjy  is  used  in  eq(2.6)  to  determine 
the  matrix  Cy.  This  matrix  in  turn  is  used  in  eq(2.3)  to  provide  the  white  noise 

driving  term  11  (n)  in  eq(2.1).  We  note  that  ii(n)  provides  the  cross-channel 
correlation  through  the  Hermitian  covariance  matrix 

In  [2],  we  describe  a  method  to  incorporate  the  desired  correlation 

properties  into  [R^].  This  is  accomplished  through  the  use  of  'shaping  functions' 

which  enable  us  to  modify  the  shape  of  the  correlation  functions  contained  in  the 
correlation  matrix.  For  the  cross -correlation  function,  we  consider 
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(2.21a) 


(pj:)a::a;if(A.:j,l  -  I::) 

R  (I)  — 'j — ^exp(j[ei/i)-e,(0)i) 


_  ipiiiqii°iif(^ii’i  -  ijj) 
-  •ij)l|=o 


(2.21b) 


where  Pjj  is  the  complex  cross-channel  correlation  parameter  such  that 
p..=R.^(0)/Ojiajj,  is  the  one-lag  temporal  cross-correlation  parameter,  ly  is  the 
lag  value  at  which  the  cross-correlation  function  peaks  and  0jj(l)  is  the  phase  of 

the  cross-correlation  function.  By  definition, 

p-  IPijI  expfj0ij(O)]  (2.22) 

Expressions  for  Rjj(l)  are  obtained  from  eq(2.21b)  using  the  property 

Rji(l)=R*(-l).  (2.23) 


These  equations  provide  us  with  a  useful  description  of  the  cross¬ 
correlation  function  in  terms  of  the  complex  cross-correlation  coefficient  Pjj,  the 
standard  deviations  Oji  and  of  the  channel  i  and  j  processes,  respectively,  and 
the  one-lag  temporal  cross -correlation  parameter,  A.jj.  For  the  autocorrelation 
function  (i=j),  we  have  lpiil=l  since  any  given  channel  process  is  totally 
correlated  with  itself  at  zero  lag.  Also,  0ii(O)  =  0  since  0ii(l)  is  an  odd  function  of 

I.  Since  the  function  f(*)  for  the  autocorrelation  function  has  a  peak  value  of  unity 
at  1=0  and  noting  that  \^^  =  0,  eq(2.21b)  for  i=j  reduces  to 

Rii(l)  =  a?f(2.ii,l)exp(jeii(l) )  (2.24) 

where  Xjj  is  the  one-lag  temporal  autocorrelation  parameter  and  is  a  measure  of 
the  correlation  magnitude  between  consecutive  samples  such  that  ()<A,jj<l.  At  lag 

value  1=0,  eq.(2.24)  becomes  Rii(O)  =  ajj  which  is,  as  expected,  the  variance  of  the 
zero  mean,  channel  i  process. 
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EXAMPLE  Gaussian  Shaped  Correlation  Functions 

In  [2],  the  special  cases  of  Gaussian,  exponential  and  sine  shaped  correlation 
functions  are  considered  using  the  forms  given  in  eqs(2.21b)  and  (2.24).  For 
cross-correlation  functions  with  Gaussian  shaped  magnitudes,  for  example,  we 


I  -  lij) = 


(2.25) 


SO  that  eq(2.21b)  becomes 


Rij(J)  =  — -  ^  fl - ^expUe,/!)} 


(2.26) 


(^ij)  'j 


For  the  autocorrelation  function  (i=j)  with  Gaussian  shaped  magnitudes,  we 
have  (dropping  the  subscript  i  for  notational  convenience) 


R(l)  =  f(>.,l)exp(je(l)]  =  (X)'  expue(l)] 


(2.27) 


where 


2 

f(X,\)  =  (k)^  =  exp[-  27c2p2x2|2] 

X  =  exp[-27c2^2q’2j 


(2.28a) 


(2.28b) 


is  a  real  constant  such  that  0  <  X  <  1  and  T  is  the  sample  period.  In  ref  [2],  we 
show  that  p2  is  the  variance  (which  determines  the  width)  of  the  corresponding 
Gaussian  shaped  spectra. 

EXAMPLE  Exponentially  Shaped  Correlation  Functions 


For  the  exponentially  shaped  cross-channel  correlation  function,  we  use  the 
expression 
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(2.29) 


Rii(l)  = 


^  /■||  ■'l - exp  j0  (I) 


while  the  autocorrelation  function  becomes 


Rii(k)  =  cfi  expljejid)]. 


(2.30) 


In  the  above  discussion,  we  have  proposed  the  use  of  functional  forms  to 
characterize  the  magnitude  and  phase  of  the  correlation  functions.  The  motivation 
for  this  approach  is  that  it  will  allow  for  flexibility  in  modeling  random  processes 
with  various  correlation  and  spectral  shape.  We  note,  however  that  at  this  point 
we  have  not  constrained  these  functions  to  meet  all  the  criteria  that  are  necessary 
and  sufficient  to  characterize  correlation  functions.  In  fact,  determining  all  of 
these  conditions  in  a  general  formulation  is  a  difficult  task.  In  [2],  we  note  several 
constraints  for  the  correlation  functions  including  the  important  condition  of 
positive  semi -definiteness  of  the  correlation  matrix. 
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3.0  VARIANCE  OF  THE  TIME-AVERAGED  CORRELATION 
FUNCTIONS 


3.1  VARIANCE  OF  THE  AUTOCORRELATION  FUNCTION  ESTIMATOR 


The  ensemble  correlation  function  is  defined  as  the  expectation  of  lagged 
products  of  a  given  stationary  process  when  averaged  over  an  ensemble  of 
realizations.  If  this  function  is  equal  to  the  time-averaged  correlation  function 
obtained  from  a  single  realization,  the  process  is  called  auto-  or  cross-correlation 
ergodic.  Consider  the  time-averaged  estimate  of  the  biasedt  correlation  function 
using  Nj  observation  time  samples 


r 


Nt-I-1 

^  n=0 


X  Xi(n)Xj(n-lll) 
n=0 


0  <  I  <N.j,-l 


-(N.J.-1)  <  I  <  0. 


(3.1) 


The  variance  of  the  autocorrelation  (j=i)  estimator  (l,Nj)  is  expressed  as 


Vb,(I,Nt)  = 

=  E  {  (I,Nt.)  -  E[feii_  (I,Nt)  (I,Nt-)  -  }  (3.2a) 

^  D  u  b  b 

=  E[^ii  (l,NT)6ii^  (I,Nt)]  -  E[Aii  (I.NT)]E[Aii*  (I.Nt-)].  (3.2b) 

‘  b  b  *  b  b 


We  now  consider  each  term  within  the  expectation  operations  expressed  in 
eq(3.2b).  Using  eq(3.1)  in  the  first  term  on  the  RHS  of  eq(3.2b),  we  have  for 
positive  and  negative  I 


Nt^-I-1  Nt-I-1 

^iiTJ*’NT)^iiTu(I.NT)  =  ^  X  X  Xj(n)xi(n-I)xi(p)xj(p-I) 

Nt  n=o  p=0 


^  In  this  paper,  wc  consider  the  biased  estimator  for  the  correlation  functions  since  it  ensures  positive  semi-definiteness.  In 
[2],  the  unbiased  estimator  is  presented  as  well  as  an  estimator  with  unlimited  data. 
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0<I<Nt-1  (3.3a) 

Nt  n=0  p=0 

-(Nt-I)^I^O  (3.3b) 

so  that 

E[ftii_,(l,NT)ftiiT,(l,NT.)]=-^  S  S  E[Xi(n)xi(n-l)x*(p)xi(p-l)] 

^  Nt  n=o  p=0 

0<I<Nt-1  (3.4a) 

Nt-III-I  Nt-II|-1  ^ 

=  -"2  S  X  E[xi(n)xj(n-lll)Xj(p)xi(p-III)] 

Nj  n=0  p=0 

-(Nt.-1)<I<0  (3.4b) 

For  the  second  term  in  eq(3.2b),  eqs.(3.1)  enable  us  to  obtain 


1 

E[^ii^(I.NT)]  =  X  Rii(l)  0  <  I  <  Nt  - 1  (3.5a) 

N>|-1, 

=  fT  X  M)  -(Nt-1)<IS0.  (3.5b) 

n=0 

SO  that 

E[ftii.jl^(l,NT.)lE[feuTb(I.NT)l  = 

Nx-l-1  Nt-I-I 

=  -^  X  X  iRiid)!^  0  <  I  <  Nt  (3.6a) 

n=0  p=0 

=  ^  I  I  •Rii(0l^  -(Nt-1)  ^1^0-  (3-6b) 

n=0  p=0 

Using  eqs(3.4)  and  (3.6)  in  the  expression 


Vb..(I,Nt.)  =  E[ftiiT5(I.NT)feiiTb(I.NT)  ]  -  E[ftiiT^(l,NT.))E[ftiiT|^(l,NT.)]  (3.7) 
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we  obtain 


Nt-I-1  Nt-I-I 

VbXI.Nt)  =  -^  I  X  (E[xi(n)xi(n-l)x*(p)Xi(p-l)]-lRii(l)|2) 

"  Nt  n=0  p=0 

for  0  <  I  <  Nt  -1  (3.8a) 

Nrll|-1NT-Il|-1  . 

=  ^  I  X  (E[xi(n)Xi(n-lll)Xi(p)Xi(p-lll)]-lRii(l)P) 

Nj  n=0  p=0 

for  -(Nt-1)  <1  <0.  (3.8b) 

We  now  define 

(t)(n,l)  =  Xi(n)xi  (n  - 1)  (3.9a) 

and 

R^^(k,l)  =  E[(|)(n,l)(t)*(n  -  k.l)]  (3.9b) 

so  that,  assuming  stationarity,  the  covariance  of  i|>(n,l)  can  be  expressed  as 

C,|^(k,l)  =  E[  { It)(n,l)  -  E[<|>(n,l)] }  { il)*(n-k,l)  -  E[(>*(n-k.l)] )  ]  (3.1  Oa) 

=  R,|„,,(k,l)  -  E[iti(n,l)lE[(|i*(n-k,l)).  (3.10b) 

Also,  from  eq(3.9a) 


E[(j)(n,l)]  =  Rii(l)  (3.11a) 

and 

E[(I)*(n-k,l)]  =  Ru(l)  (3.11b) 

so  that  (3.10b)  becomes 

C,,^(k,l)  =  R,f^(k,l)  -  lRii(l)P  (3. 1 2a) 

=  E[xj(n)xj(n  -  l)xi(n  -  k)xj(n  - 1  -  k)|  -  lRji(l)P.  (3.12b) 


Using  eq(3.12b)  in  (3.8) 
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(3.13a) 


Nt-I-1  N-r-l-l 

Vb,(^n^)  =  ;:;2  2  2  c^4>(n  -  p,i)  o  <  i  <  - 1 

n=0  p=0 

N-t-III-INt^-III-I^ 

=  ^2  S  I  C..(n  -  p,lll)  -(Nt-1)<I<0.  (3.13b) 

Nt  n=0  p=0 

Using  k  =  n  -  p  where 


(Nt- - 1  -  1)  <  k  <  - 1  -  1 

for 

0<l  <Nt-1 

(3.14a) 

(Nt--II1-  l)<k<NT,-lll-  1 

for 

-(Nt.-1)<1<0. 

(3.14b) 

We  also  note  that  eq(3.14b)  holds  for  all  I  so  that 


Nt-I-1 

VbhC.Nt)  =  ^2  I  [Nt  -  III  -  Ikl]  C.*(k,l)  0  <  I  <  Nt  - 1  (3. 1 5a) 

"  Nx  k=.(Nx.|-l) 

Nx-ll|-l 

=  ;:!2  X  [Nx  -  III  -  Ikl]  cL(k,lll)  -(Nx-l)<l<0.  (3.15b) 

Ny  k=-(Nx-ll|-l) 

However,  for  negative  lag  I,  we  have 


C„,^(k,lll)  =  C^(k,l)  (3.16) 

and  eqs(3.15)  can  be  written  as 


Vb.C.Nt)  = 


Tk=-(NT.-ll|-l) 


III  +  Ikl 
Nx 


C(])(j,(k,l) 


(3.17) 


for  both  positive  and  negative  values  of  I.  For  processes  with  zero-mean,  jointly 
stationary  Gaussian  quadrature  components,  the  imaginary  terms  in  C(j)(])(k,l) 

cancel  when  summed  over  positive  and  negative  values  of  k  (Appendix  A)  so  that 
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1  r  iii+iki 

Vb  (I,Nt)=  jq-  I  1  -  Re{C^(k,l)l. 


(3.18a) 


if  S  1  [lRij(k)P  +  Re(Fii(l,k))]  (3.18b) 

T  k=-(NT,-ll|-l)'-  •' 


where 


Rij(k)  =  E[Xi(n)xi(n  -  k)] 


(3.19) 


Fji(l,k)  =  E[Xj(n)Xi(n  -  I  -  k)]E[xi(n  -  l)xi(n  -  k)].  (3.20) 

Using  eq(A.6b)  in  Appendix  A,  it  can  be  shown  that  in  the  special  case  where  the 
quadrature  components  maintain  the  relations 


Ru(l)  =  R^(l) 


(3.21a) 


r!P(I)  =  -  Ra‘(l) 


(3.21b) 


Re(Fii(l,k))  =0 


for  all  l,k. 


(3.22) 


We  note  that  eqs,(21)  hold  for  the  special  case  of  wide-sense  stationarity  of  the 
narrowband  bandpass  processes  [2].  In  this  case,  eq(3.18b)  becomes 


N^-lll-1 


VBi-O.N^)  =  r  S  [l  - 

'rk=.(N  ,111-1)1-  T’-l 


(3.23) 


3.2  VARIANCE  OF  THE  CROSS-CORRELATION  FUNCTION  ESTIMATOR 


Analogous  to  the  derivation  leading  to  eq(3.18b),  the  variance  of  the 
biased,  time-averaged  cross-correlation  function  estimator  for  processes 
described  by  the  general  Gaussian  case  noted  in  section  3.1  has  been  shown  (2]  to 
be 


Nt-III-I 


VBij(I.NT)  =  I 

^  ^k=-(NT-lll-l) 


1  - 


lll+lkl 


N. 


T  J 


Re[Rii(k)Rjj(k)  +  Fij(l,k)]  (3.24) 


where 


Fij(l,k)  =  E[Xi(n)Xj(n  - 1  -  k)]E[xi  (n  -  l)xj  (n  -  k)].  (3.25) 

We  note  that  the  term  Re[Fjj(l,k)]  will  contribute  a  dependence  of  Vg..(l,N)  upon 
cross-correlation  terms  such  as  the  cross-correlation  parameter  Ipjjl  defined  in 

section  2.4  and  [1,2].  In  [2],  we  show  that  for  joint  wide-sense  stationarity  of 
narrowband  multichannel  bandpass  processes, 

r“(I)  =  R^^d)  (3.26a) 

and 

Rif  (I)  =  -  Ry ‘(I)  (3.26b) 

which  leads  to  Fjj(l,k)=0.  Under  these  conditions,  eq(3.24)  reduces  to 


1 


VBij(I.NT)  =  N 


Nt-III-1 

I 


1  - 


Re[Rii(k)R^(k)]. 


(3.27) 


k=-(NT.-lll-l)L 

Thus,  for  joint  wide-sense  stationarity  of  the  narrowband  multichannel 
bandpass  processes,  the  variance  of  the  cross-channel  correlation  function 
estimator  is  independent  of  the  cross-channel  correlation.  This  result  is  verified 
via  simulated  results  in  sections  5.2  and  5.4. 


I 


3.3  SAMPLE  VARIANCE  OF  THE  CORRELATION  FUNCTION 

Consider  Nj^  realizations  of  the  random  process  Xi(n).  Let  each  realization 
be  indexed  by  the  integer  a;  a=l,2,...,Nj^.  Corresponding  to  the  realization  with 

index  a,  let  ]^ijjj^(l,Nj|a)  be  the  biased,  time-averaged  cross-correlation  function 
estimate  using  Nj  observation  samples.  The  sample  variance  of  the  time -averaged 
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cross-correlation  function  estimate  is  computed  from  Nj^  statistically  independent 
realizations  using  the  expression 


1  — 

Var[ftijTj,(I.NT):NR]  =  N-lj  S  I  iJ  ijr.  (l,NT|a)  -  fijjT-  (l,NRla)P  (3 .28) 

^  a=l 

where 

Nr 

ftijT^(I.NRla)  =  I  ftijTjl.NTl«)-  (3.29) 

*^a=l 

The  autocorrelation  function  is  obtained  for  i=j.  Eq.(3.28)  is  used  to  compute  the 
variance  of  the  time-averaged  auto-  and  cross-correlation  function  estimators  in 
chapter  5  for  comparison  with  eqs(3.23)  and  (3.27),  respectively. 


3.4  ERGODICITY  CONSIDERATIONS 


In  section  2.4  and  [1,2],  we  presented  functional  forms  for  the  auto-  and 
cross-correlation  functions  to  obtain  desired  temporal  and  cross-channel 
correlation  of  the  processes  to  be  synthesized.  Using  eq(2.30)  in  eqs(3.23)  and 
(3.27),  we  obtain  [2],  respectively 


1 


Nt-III-I 


1  - 


lll-Hlkl 
Nt  j 


Oii  (>ni) 


21kl 


(3.30) 


and 


VBiid.N^)  = 


Nt^-III-I 


N. 


T  V- 


1  - 


k=-(NT-lll-l)'- 


lll-flkl 
N. 


T  J 


J-  n  \lkl  2  .y  Jkl  . 
^11  V^n)  ^22  (^22^ 


cos[eij(i)  -  ejj(i)]. 


(3.31) 


When  invoking  the  ergodic  assumption,  one  would  like  the  variances  of  the 
estimators  expressed  by  eqs.(3.30)  and  (3.31)  to  be  suitably  small  so  that  the 
time-averaged  correlation  function  is  a  satisfactory  approximation  to  the 
ensemble  correlation  function.  For  a  given  Ajj  and  these  analytic  expressions 
provide  a  means  for  determining  the  required  Nj  to  minimize  the  variance  to  a 

specified  level.  A  unique  aspect  of  this  development  is  the  determination  of  these 
expressions  in  terms  of  both  the  observation  window  size  Nj,  the  process 

parameters  ajj  and  and,  for  V3.j(l,Nj),  the  phase  angles  0jj(l)  and  0jj(l).  These 
expressions  indicate  that  for  ^ii<l,  their  limit  approaches  zero  as  Nj  approaches 
infinity.  Thus,  for  stationary  processes,  ergodicity  holds  in  all  cases  except  for 
total  temporal  correlation  (ie..  A,;—!).  In  this  case,  Vg  (l,Nj)=Ojj;  ie.,  the  square 

of  the  process  variance.  Furthermore,  as  of-  and  Ajj  change,  eqs(3.30)  and  (3.31) 

provide  a  quantitative  measure  expressing  the  requirements  on  the  observation 
window  size  to  obtain  time-averaged  correlation  function  estimates  which  yield  a 

close  approximation  to  the  ensemble  averaged  correlation  functions. 

Figure  1  is  a  plot  of  the  error  variance  Vg  .(l,N-p)  for  the  biased,  time- 

averaged  autocorrelation  function  estimator  versus  the  number  of  time  samples 
Nj  used  in  the  estimate.  Each  curve  is  plotted  for  a  specific  value  of  the  temporal 
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2 

correlation  parameter  X-  while  Oji  was  fixed  at  a  value  of  4.  As  these  curves 
reveal,  the  error  variance  decreases  with  increasing  N-p,  however,  the  rate  of  the 

decrease  is  highly  dependent  upon  the  temporal  correlation  of  the  processes;  ie., 
for  processes  with  low  temporal  correlation  (^jj  approaching  zero),  V3,  (l,Ny) 

diminishes  rapidly  for  increasing  N-p.  However,  for  processes  with  high  temporal 
correlation  (Xjj  approaching  unity),  Vg..(l,Ny)  decreases  slowly.  At  1,  the 

process  is  no  longer  ergodic,  so  that  the  error  variance  no  longer  decreases  as  a 
function  of  N-p. 


>.=0.5 

>.=0.9 

>.=0.99 

>.=0.999 

h=0.9999 

>1=1.0 


Figure  1  Variance  Vg..(l,Nj)  at  1=0  for  the  time-averaged  autocorrelation 

function  estimator  as  a  function  of  Nj  time  samples  with  the  one-lag 

2 

temporal  correlation  >.„  as  a  parameter  and  aii=4. 
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3.5  THE  NORMALIZED  CORRELATION  FUNCTION 


The  normalized  time-averaged  auto-  and  cross-correlation  function 
estimators  are  briefly  discribed  in  this  section.  They  are  defined,  respectively,  as 

fii(l,NT)  =  ftii(I.NT)/ftii(0,NT)  (3.32) 

and 

(\j(l,NT)  =  ftij(I.NT)/ftij(0,NT.).  (3.33) 

Specific  forms  of  the  estimators  such  as  the  biased  correlation  function  estimator 
presented  in  eq(3.1)  are  used  to  evaluate  the  numerator  and  denominators  of 
eqs.(3.32)  and  (3.33).  The  normalized  forms  differ  significantly  from  the 
corresponding  unnormalized  form  presented  above  in  as  much  as  the  variance  of 
the  estimates  have  a  significantly  different  dependence  upon  process  correlation. 
We  note  that  each  estimator  is  expressed  as  a  ratio  of  two  random  variables. 
Each  of  these  random  variables,  however,  is  highly  correlated;  ie.,  they  are  each 
computed  using  the  same  observation  data.  As  a  result,  the  deviation  of  ftij(l,Nj) 

about  its  mean  level  will  be  'tracked'  similarly  by  ftij(0,Nj);  ie.,  if  ftij(l,Ny)  lies 

above  (or  below)  its  mean  level,  likewise,  ^jj(0,NT)  will  lie  above  (or  below)  its 

mean  level.  As  a  result,  the  normalized  correlation  function  estimator  will  be 
closer  to  its  respective  mean  and  its  associated  variance  is  reduced.  As  noted 
above,  its  dependence  upon  process  correlation  will  be  quite  unlike  that  of  the 
ordinary  correlation  function  estimators.  In  fact,  it  will  be  similar  to  that  of  the 
parameter  estimators  described  later  (see  section  4.5). 

The  variances  associated  with  the  estimators  described  by  eqs(3.32)  and 
(3.33)  are  difficult  to  develop  since  they  involve  ratios  of  two  correlated  random 
variables.  In  [17],  a  taylor  series  expansion  is  used  to  derive  approximate 
expressions  for  the  variance  of  these  estimators. 
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4.0  AUTOREGRESSIVE  PROCESS  PARAMETER  ESTIMATORS 


In  this  chapter,  we  describe  the  three  parameter  estimators  used  in  this 
study;  namely,  the  multichannel  Strand-Nuttall  estimator  [11,12]  and  two  versions 
of  the  Yule-Walker  estimator  [9,10].  In  addition,  we  present  a  recursive  method 
used  to  calculate  the  multichannel  matrix  coefficients  known  as  the  Levinson- 
Wiggins-Robinson  method.  The  Burg  estimator  [13]  and  the  multichannel 
Strand-Nuttall  estimator  are  equivalent  for  the  single  channel  case. 

Expressions  for  the  error  variance  of  time-averaged  parameter  estimators 
are  difficult  to  determine.  However,  the  Cramer-Rao  bound  (CRB)  provides  a 
reference  for  the  performance  of  the  estimators.  An  analytic  expression  for  the 
exact  Cramer-Rao  bound  (CRB)  of  unbiased  AR  parameter  estimates  has  been 
presented  [6].  The  expression  was  shown  to  converge  to  the  asymptotic  form  of 
the  CRB  for  large  measurement  time  sample  sizes.  In  chapter  5.0,  we  show 
simulated  results  for  the  error  variances  of  the  estimators  and  compare  these 
results  to  the  exact  CRB  as  a  function  of  the  correlation  of  the  processes. 

4.1  THE  YULE- WALKER  ESTIMATOR 

The  Yule-Walker  estimator  used  in  this  study  initially  uses  estimates  in  the 
covariance  matrix  of  eq.(2.14).  The  Yule-Walker  equation  is  then  solved  to 
obtain  the  estimates  of  the  parameters  and  the  white  driving  noise  covariance 
matrix.  The  two  versions  of  the  Yule-Walker  estimators  used  in  this  study  differ 
in  that  the  estimates  of  the  covariance  matrix  are  obtained  using  either  the  biased 
or  the  unbiased  covariance  estimator.  They  are  refered  to  here  as  the  Yule- 
Walker  with  biased  covariance  (YWBC)  estimator  and  the  Yule-Walker  with 
unbiased  covariance  (YWUBC)  estimator,  respectively.  The  biased  covariance 

estimator  is  expressed  in  eq.(3.1).  The  corresponding  unbiased  estimator  is 
obtained  by  replacing  the  denominator  of  eq.(3.1)  by  Nj-I.  We  point  out, 

however,  that  although  the  covariance  estimator  may  be  unbiased,  the 
corresponding  Yule- Walker  (YWUBC)  estimator  is  not  unbiased. 

4.2  THE  LEVINSON- WIGGINS-ROBINSON  RECURSION 

We  now  develop  a  recursive  procedure  to  determine  the  multichannel 
forward  and  backward  matrix  coefficients.  The  procedure  was  initially  presented 
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by  Wiggins  and  Robinson  [16]  based  on  the  scalar  version  of  Levinson  [14].  The 
approach  presented  in  this  subsection,  however,  is  similar  to  that  used  by  Strand 
[11].  Consider  the  pth  order  update  of  the  multichannel  forward  PEF  coefficient 
expressed  in  the  form 


II 

o. 

<1 

"  Ap(0)  - 
Ap(l) 
Ap(2) 

■  Ap.,(0) 
Ap-,(1) 
Ap-i(2) 

+ 

[0] 

Bp-i(p-l) 

Bp-,(2) 

Ap(p-l) 

-  Ap(p)  _ 

Ap-,(P-1) 

-  [0]  - 

Bp-,(1) 

_  Bp.,(0)  _ 

(4.1) 


where  [Fflp  is  a  constant  matrix  called  the  forward  reflection  coefficient  matrix 
and  Ap(0)=Ap.j(0)=Bp.i(0)  =  I,  the  JxJ  identity  matrix.  TTie  update  equation  is 
therefore  expressed  as 


=  Ap.,(k)  +  Bp.,(p-k)[rf]p 

k=l  ,2,...,p-l 

(4.2a) 

=  [rflp 

k=p 

(4.2b) 

=  [0] 

k>p 

(4.2c) 

=  Ap(0)  =  I 

k=0. 

(4.2d) 

Consider  the  order  update  of  the  backward  PEF  coefficient  as 

~  Bp(P)  ~ 

Bp(p-l) 

Bp=  : 

Bp(l) 

_  Bp(0)  _ 


[0] 

Bp-, (P-1) 

Ap.,(0) 

Ap.,(l) 

Bp-,(1) 

_  Bn.,(0)  _ 

+ 

Ap-,(P-1) 

_  101  _ 

(4.3) 


where  [rfj]p  is  the  backward  reflection  coefficient  matrix.  The  update  equation  is 
expressed  as 


Bp(k)  =  Bp.,(k)  +  Ap.,(p-k)[rblp  k=l,2,...,p-l 


(4.4a) 
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=  [^blp 

k=p 

(4.4b) 

=  [0] 

lop 

(4.4c) 

=  Bp(0)  =  I 

k=0. 

(4.4d) 

We  now  determine  the  conditions  under  which  the  above  equations  hold. 
First,  we  take  the  Hermitian  transpose  of  eq(4.1)  and  postmultiply  both  sides  of 

the  equation  by  [Rxx]p+i-  LHS,  the  normal  equations  provide 

[Rm'p+i  =  ( [“1  •  [01  }  (''•5a) 

PJ 

where  [Lf]p  is  the  JxJ  forward  error  covariance  matrix  for  a  pth  order  filter. 

Each  side  of  eq(4.5a)  is  a  vector  of  matrices  containing  p+1  matrix  elements. 
Taking  the  Hermitian  transpose  of  Ap  in  eq(4.1)  and  using  it  in  (4.5a),  we  also 

have 

Ap  [Rxxlp*.  =  {  aJI,  [0]  }  +  [Ffl^  {  [0]  }  (Rjj^lp,,  (4.5b) 

SO  that 

{ [Iflp  [0] ...  [0|  1  =  {  A^.,  [0] }  (Rjyjlp^i  +  [fflp  {  [0]  Bp., }  [Rxxlp.,.  (4.5c) 
Let  us  now  introduce  a  partitioned  form  of  [Rxxlp+i 

R„(0)  R„(l)  R„(2)  ...  R„(p) 

R„(-1)  R.x(0)  Rxx(i)  •Rxx(P-i) 

•  •  »  •  • 

•  •  •  •  • 

Rxx(l-P)Rxx(2-P)Rxx<3-P)  --  Rxx(l) 

Rxx(-P)  Rxx(i-P)Rxx(2-P)  -  RxxO) 

where  (Sxx)p  is  a  vector  containing  p  matrix  elements  and  is  defined  as 

(Sxx)?  =  [Rxx(-P).Rxx(i-P).-.Rxx(-i)l-  (4-7) 

We  now  consider 


[Rxxip  ^xx^p 
-  P  ^  (4.6) 

(Sxx)plRxx«))lJ 
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(Ap-i  [0]  )  (Rjyilp,.!  -  {Ap-1  [01  1  (4.8a) 

L(&xx)^Rxx(0)l. 

=  {a",[R^1p  A^.,(Sxx)p}-  (4-8b) 

However,  we  can  reduce  the  order  of  eq(4.5a)  by  one  and  obtain  the  Hermitian 
transposed  vector  which  will  now  contain  p  matrix  elements,  so  that 

A^iPjtxlp  =  (IJ^lp-i  [01  -  [Oil.  (4.9) 

Let  us  now  define 

[Al"  =  A^.i(Sxx)p  =  [I  Ap.,(l)  Ap,|(2) ...  Ap.,(p-1))  r  Rxx(P)  1  (4.10a) 

Rx,(P-i) 

L  R„(l)  _ 

=  Rxx(P)  +  L  Ap.,(k)R„(p-k).  (4.10b) 

k=l 

Substituting  eqs(4.9)  and  (4.10a)  in  (4.8b),  we  have 

[A^-i  [Oil  [R^j^lp+i  =  [[Ilflp.,  [01 ...  [01  [Alp  1.  (4.11) 


We  will  utilize  this  expression  later.  We  now  repartitioned  [R^]p+i  into 
the  form  given  by 

■Rxx(0)  Ml)  R,x(2)  Rxx(p)  ■ 

Rxx(-i)  RxxW)  RxxC)  •  ■•Rxx(P-')  rRj^(0)(Rxx)p' 

[Rjilp+i  “  •  •  :  :  :  “  .  h?  i  (4.12) 

R„(l-p)R„(2-p)R„(3-p)...  R„(l)  L(Rxx)p  [Rxxlp-I 

-Rxx(-P)Rxx(l-P)Rxx<2-p)-  Rxx(O)  - 

where 

(Rxx)p  =  (Rxxd)  Rxx(2)  •••  R,x(p)]-  (4.13) 

We  now  consider 
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{[0]  }(R^lj,^,  =  {[01  Bp-i}rR,,(0)(R,,)"l  (4,14a) 

-  (^xx)p  P^XX^p  - 

=  {B^.,(R„)p  a”,[Rii)pl-  (414b) 

However,  from  the  backward  multichannel  normal  equations,  we  have 


4^i[RMlp=f(01-  [01  l^blp-il- 


Let  us  now  define 


lA'ip  =  a^.,  (R,,)p 

=  {Bp.,(p-l)Bp.|(p-2) 


•  Bp., (1)1)  rRxx(-i) 

Rxx(-2) 


Bxx(-P) 


=  R„(-p)  -b  £  Bp.,(k)R,j(k-p). 

k=l 


Substituting  eqs(4.15)  and  (4.16a)  into  (4.14b),  yields 


(4.15) 


(4.16a) 

(4.16b) 


(4.16c) 


{ (01  b“,  )  (Rj^]p^,  =  { (A’l"  [0] ...  [0]  [EbJp-i )  • 

Using  eqs(4.1 1)  and  (4.17)  in  (4.5c),  we  obtain 


(4.17) 


([Sf]p  [01 ...  (Oj)  =  {l£fjp.|  [0] ...  [0]  [Alp  ) 


p-1  •••  iv/j  ItJJp  J 

+  [rfip  {(A'lp  [01  ...[01  dbip.,}. 


+  [Pflp  1  [A'lp  [01 ...  [01  [Sblp., ) .  (4.18) 

And  so,  in  order  for  eq(4.1)  to  hold,  we  must  satisfy  eq(4.18).  From  eq(4.18), 
we  have 


[^f]p  “  [^f]p-l  [^f]p  ]p 


(4.19) 


[01  =  [Alp  +  [Pflp  [Zblp-I- 

From  eqs(4.20)  and  (4.2b),  we  have 


(4.20) 
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[Hflp  =  Ap (p)  =  -  [Alp  { Hblp-i )■'•  (4-21) 

We  must  now  return  to  eq(4.3)  and  carry  out  a  similar  procedure  that  was 

H 

carried  out  for  eq(4.1);  ie.,  we  post-multiply  Bp  by  [Rxx]p+i  that 

bJ*  [Rxxlp+i  =  f  [01  [01  l^blp  }  ■  (4.22) 

For  the  first  term  on  the  RHS  of  eq(4.3),  we  postmultiply  by  the 
partitioned  form  of  eq(4,12).  In  this  case,  we  obtain  the  same  results  as 
expressed  in  eqs(4.14a)  through  (4.17).  For  the  second  term  on  the  RHS,  we 
postmultiply  by  the  partitioned  form  of  eq(4.6)  so  that  we  obtain  the  same  results 
as  expressed  in  eqs(4.8)  through  (4.11).  Substituting  these  results  into  eq(4.22), 
we  have 


26 


{ [01 ...  [0]  [£blp  1  =  { [A')p  [0| ...  [0]  Hblp-I  ) 

+  [Fblp  {[Eflp.,  [01 ...  [01  [Alp  1 .  (4.23) 

Eq(4.23)  must  therefore  be  maintained  in  order  for  the  recursion  in  eq(4.3)  to 
hold.  From  the  above  equation,  we  have 

[Zblp  =  [Sblp.,  +  [fblp  [Al" 

[01  =  [A’lp  +  [Fblp  [Xflp-,. 

From  eqs(4.25)  and  (4.4b),  we  have 

irblp  =  Bp  (p)  =  -  [A'lp  {  [Sfip.,  )■'.  (4.26) 

Before  summarizing  the  pertinent  recursion  equations  used  to  estimate  the 
forward  and  backward  coefficient  matrices,  we  can  at  this  point  develop  several 
important  relationships.  Using  equations  (4.21)  and  (4.25)  in  (4.19),  we  obtain 


[Iflp  =  [Zflp.,  +  [Pflp  [A'lp 

(4.27a) 

=  [Sflp,,  -  a"(p)  [fblp  [Sflp.,. 

(4.27b) 

Using  eq(4.26)  in  (4.27b),  we  have 

,H  ,H  H 

[^f]p  ~  [^f]p-l  ■  ^p  (p)  ®p  (p)  [^f]p-l 

(4.28a) 

so  that 

[2;fl“={l-Ap(p)Bp(p)}  [Zflp.,. 

(4.28b) 

Similarly, 

[Zblp  =  {l-Bp(p)Ap(p))  [Ibip-I. 

(4.29) 

(4.24) 

(4.25) 


Let  us  also  consider  the  matrix  identity 
,  .  H 


1011  IRxxlp+i 

[or 

[fip-.J 

H 


{(01  Bp.,}  [R^jp^, 


^p-i 

L[0]  J 


H 


(4.30) 
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which  follows  since  [Rxxlp+i~  t^xxlp+i*  Using  eqs(4.11)  and  (4.17)  in  (4.30),  we 
can  write 


{[5:f]p.l[0]...[0][A]»}p  foj  - 

Bp-l(p-l) 

Bp.i(l) 


L  I 


•"{[A’lp  [0]  .. 


I0][2:b]p.i)p  I 

Ap.i(2) 

Ap.i(p-l) 

L  [0]  -1 


H 


(4.31) 


so  that  after  carrying  out  the  matrix  multiplications  in  eq(4.31),  we  have 

[Alp  =  [A'lp.  (4.32) 

Eq(4.32)  is  useful  since  it  reduces  the  computations  involved  in  the 
determination  of  the  reflection  coefficient  matrices;  ie.,  eqs(4.21)  and  (4.26)  can 
both  be  expressed  in  terms  of  [A]p,  so  that 

[Fflp  =  Ap(p)  =  -  [Alp  {[41p-i}''  (4.33a) 

[Ftlp  =  Bp  (p)  =  -  (Alp { [Zfl" , )  - ‘ .  (4.33b) 

In  the  single  channel  case,  the  scalar  forward  and  backward  prediction 

error  variances  are  equal.  Since  [A]p  is  a  scalar,  it  follows  from  eqs(4.33)  that 
* 

ap(p)=bp(p). 

We  again  note  that  the  expressions  developed  here  are  in  terms  of  the 
Hermitian  operator  as  opposed  to  those  often  shown  in  the  literature.  This  stems 
from  the  notation  used  to  express  the  autoregressive  equation  in  eq(2.1)  and 
provides  a  more  consistent  notation  as  stated  previously.  We  now  summarize  the 
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pertinent  equations  for  the  recursion  to  update  the  forward  and  backward  PEF 
coefficients.  First,  we  have  the  definitions 


[A]p  =  +  X  Ap.,(k)R,,(p-k)  (4.34a) 

k=l 

[A']p  =  R„(-p)  +  X  Bp.,(k)R„(k-p).  (4.34b) 

k=l 

However,  as  shown  previously 

[A]"  =  [A']p.  (4.35) 

Hence,  we  need  use  only  one  of  the  equations  in  eqs(4.34).  Next,  from  eqs(4.33) 


irrfp  =Ap(p)  =  -[AJp{[Eblp.,)‘' 

[Fblp  =Bp(p)  =  -[Alp{[Ef]p.,)-' 

where 

[2f)p  =  ( I  -  Ap  (p)  B^(p)} 

(4)p  =  {l-Bp(p)Ap(p)}  Pblp.,. 


Then, 

Ap  (k)  = 
Bp  (k)  = 


'Ap.,(k)+(rflpBp.,(p-k)  k=  1 , 2,  ...,  p- 1 

< 

irflp  k=p 

'Bp.,(k)  +  [fblp  Ap.,(p-k)  k=l,  2 . p-1 

< 

[r^ip  k=p. 


(4.36a) 

(4.36b) 

(4.37a) 

(4.37b) 

(4.38a) 

(4.38b) 


We  initiate  the  procedure  using  the  initial  conditions  for  eqs(4.37)  as 
=  Rxx(O)  so  that  eqs(4.38)  provide 


A“(l)  =  [rf]“  =  -[All‘  (lSbl“)'' 

(4.39a) 

=  -R„(l)R,',(0) 

(4.39b) 

Br(l)  =  irb)r  =  -[Ali  (Wo”}  ' 

(4.39c) 
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(4.39d) 


=  -R^^(-l)Ri(0). 

The  procedure  continues  recursively  until  the  desired  order  is  reached. 

4.3  THE  STRAND-NUTTALL  METHOD 

The  multichannel  harmonic  mean  algorithm  which  is  a  generalization  of 
the  Burg  method  was  proposed  independently  by  Strand  [11]  and  Nuttall  [12].  In 
the  single  channel  case,  Burg  [13]  estimated  the  reflection  coefficients  directly 
from  the  data  subject  to  the  constraint  of  minimizing  a  performance  index 
consisting  of  the  mean-square  values  of  the  estimates  of  the  forward  and 
backward  prediction  errors.  In  this  case,  the  scalar  forward  coefficients  of  linear 
prediction  are  the  complex  conjugates  of  the  backward  coefficients.  In  the 
multichannel  case,  however,  the  forward  and  backward  reflection  coefficient 
matrices  are  related  by  eqs(4.33a)  and  (4.33b),  respectively.  These  equations 
indicate  that  the  multichannel  forward  and  backward  reflection  coefficient 
matrices  have  a  more  complicated  relationship  than  in  the  scalar  case. 
Specifically,  they  are  related  through  the  forward  and  backward  error  covariance 
matrices.  Solving  eqs(4.33)  for  -[A]p  and  noting  that  the  error  covariance 
matrices  are  Hermitian,  we  have 

-[Alp  =  Ap(p)  [Sblp.,  =  [Sflp.|  Bp(p).  (4.40) 

From  eq(4.38)  at  k=p,  we  have 

-  [Alp  =  irflpUbl",  =  [Zflp-I  [Pblp.  (4.41) 

These  equations  show  the  relationship  between  the  forward  and  backward 
reflection  coefficients  and  error  covariance  matrices  as  well  as  the  matrix  [A]p . 
In  the  Strand-Nuttall  algorithm,  it  is  this  latter  matrix  that  we  solve  for.  In 
eq(4.41),  note  that  the  error  covariances  can  be  obtained  from  stage  p-1.  Two 
additional  well  known  relationships  between  the  forward  and  backward  error 
residuals  can  be  expressed  as  [10,1 1] 
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and 


(4.42) 


6p(n)  =  £p.,(n)  +  Ap  (p)  Bp.,(n-1) 

fip(n)  =  fip.i(n-l)  +  Bp  (p)  £p.,(n).  (4.43) 

In  this  section,  we  will  establish  a  performance  function  which  depends 
upon  estimates  of  the  forward  and  backward  error  covariance  matrices.  This 
performance  function  will  be  minimized  with  respect  to  the  matrix  [A]p  .  In  the 

multichannel  case,  the  reflection  coefficients  are  chosen  such  that  eqs(4.41) 

through  (4.43)  are  satisfied.  These  are  the  constraint  equations  under  which 
-[A]p  will  be  selected  to  minimize  the  performance  function.  Using  the  notation 

of  Nuttall  [12],  we  consider  the  unbiased  error  residual  matrices  over  the 
available  data  for  a  filter  of  order  p  as 

Ep  =  S  £p(k)£p  (k)  =  Ep  (4.44a) 

^  k=p+l 

Fp  =  S  &p(k)ap  (k)  =  F^.  (4.44b) 

^  k=p+l 

Using  eqs(4.42)  and  (4.43)  in  the  above  equations  and  expanding  yields 

Ep  ^  I  [£p-,(k)£p-,(k)  +  Ap(p)fip.,(k-l)£^.,(k)  + 

^  k=p+l 

£p-i(k)fip-i(k-l)Ap(p)  +  Ap(p)fip.,(k-l)fip.,(k-l)Ap(p)]  (4.45a) 

and 

Fp ^ 5^  I  [ap.,(k-i)fl“,(k-i) + Bp(p)£p.,(k)ep-,(k-i) + 

^  k=p+l 

ep.,(k-l)E^.,(k)Bp(p)  +  Bp  (p)£p.,(k)£j‘,(k)Bp(p)]  (4.45b) 

SO  that 

Ep  =  [ftfflp-i  +  Ap  (p)  (ftfblp-i  +  (ftn>lp-,Ap(p)  +  Ap  (p)[ftbblp.,Ap(p)  (4.46a) 

Fp  =  [^bblp-l  +  Bp  (P)  (ftfb)p-i  +  [ftfbl“iBp(p)  +  Bp  (p)|ftfflp.,Bp(p)  (4.46b) 

where  the  estimated  time-averaged  error  covariance  matrices  are  expressed  as 
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(4.47a) 


[ftfflp-.  =  i  £p.,(k)£p.,(k) 

^  k=p+l 

[^fl>lp-i=[^  S£p-i(k)ap-,(k-l)  (4.47b) 

^  k=p+l 

l^fblp-i=i;^  Sap-l(k-l)£"i(k)  (4.47c) 

^  k=p+l 

[^bbip-i=s^  <‘*-4^‘*) 

^  k=p+l 

In  the  Strand-Nuttall  algorithm,  we  minimize  the  sum  of  the  traces  of  the 
weighted  error  matrices 

Isn«Alp)  =  tr(Ap.,Ep)  +  tKOp-iFp)  (4.48a) 

where  Ap.^  and  Qp.j  are  Hermitian,  positive  definite  matrices  which  weight  the 
matrices  Ep  and  Fp,  respectively.  These  matrices  provide  arbitrary  weights  and 

their  selection  will  be  determined  later.  Eq(4.48a)  is  the  performance  function 

mentioned  previously.  In  the  single  channel  case,  we  can  consider  minimizing  an 
index  of  performance  Ip(rp)  expressed  here  as  the  weighted  sum  of  the  mean- 

squared  values  of  the  forward  and  backward  prediction  errors  ep(n)  and  |3p(n), 
respectively;  ie., 

yrp)  =  aE[lEp(n)|2]  +  (I-a)E[iPp(n)l^].  (4.48b) 

Burg  considered  the  minimization  of  eq(4.48b)  for  the  case  where  a=l/2. 
We  will  show  that  the  minimization  of  eq(4.48a)  as  developed  here,  leads  to  a 
generalization  of  the  results  obtained  by  the  minimization  of  eq(4.48b)  with 
a=l/2.  Using  eq(4.40),  let  us  define  Gp  such  that 

Gp  H  -  [A]p  =  Ap  (p)  [Eblp.,  =  [Zflp.,  Bp(p).  (4.49) 

Changing  notation  for  the  Hermitian  error  covariance  matrices  to  a  form  similar 
to  Nuttall's  [10],  let 
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Up.i  =  [Ef]p.,  (4.50a) 

Vp.i  =  [Sblp.^  (4.50b) 

so  that 

Ap  (p)  =  Gp  [Vp.,]H  (4.51a) 

Bp  (p)  =  Gp[Up.,lH  (4.51b) 


Using  the  last  two  equations  together  with  eqs(4.46),  we  can  write 

Ap.,Ep  +  np.,Fp  =  Ap.,  { [ftfflp.,  +  Gp  [Vp'.,]H[ftfl,]”,  +  [ftn,lp-i[V^]Gp 

+  Gp[Vp‘.,]H[ft,,Jp.,[Vp'.,]Gp} 
iip-i  {  [^bblp-i  +  Gp[U-;,]H[ftfl,]p.,  +  [ftfl,l“,[Up'.,]Gp 

+  Gp[U',J.,]H[A„]p_,[U,i,lGp  }  .  (4.52) 

Taking  the  trace  of  eq(4.52)  and  rearranging  terms,  we  obtain 

tr[Ap.iEp  +  ^ip.^Fp]  = 

=  tr{  Ap.,[ftfflp.,  +  £2p.,[ftbblp-i  +  Ap.,Gp  [Vp'.,]H[ftn,]” 

+  £2p-i[^fbl“i(G^lGp} 

+  tr{  [Ap-i(ftfblp.,[Vp'-,]  +  Ap.,Gp  [V;;.,]H[ftb5]p.,[Vp'.,)]Gp 

+  np.,Gp[Up'.,]H(ftfl,)p.,  +  np.,Gp[Up'.i]H[ft„]p.,[Up'..lGp  }  . 

(4.53) 

We  now  consider  the  matrix  relation 

tr[KL]  =  tr[LK]  (4.54) 

where  K  and  L  are  conformable  matrices.  Noting  that  all  the  above  matrices  are 
JxJ,  the  last  two  terms  in  the  second  trace  of  eq(4.53)  can  be  reexpressed  (ie., 
with  K=  Qp.jGp  and  L  defined  accordingly  for  each  of  these  last  two  terms).  Gp 

can  then  be  factored  out  in  the  second  trace  term  so  that 
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tr[Ap-jEp  +  Qp.jFpl  - 

=  tr{  Ap.,[ftfflp.,  +  £ip.,lftbblp-i  +  Ap.,Gp  [Vp‘ 

+  tr{  [Ap.,lftfblp.,[Vp‘.,l  +  Ap.,Gp  lVp'.,]H[ft^^lp.,(Vp'.,l 

+  IUp‘.,lH|ftfblp.,Gp.i  +  IU;i,l»lAfflp.,[Up‘.,)Gp  )np.,]Gp  }  . 

(4.55) 

H 

Nuttall  now  minimizes  eq(4.55)  by  choice  of  Gp  (in  the  Hermitian  notation 
used  here),  subject  to  eqs(4.41)  through  (4.43).  Eq(4.49)  can  then  be  used  to 

H  H 

solve  for  the  coefficients  Ap(p)  and  Bp(p)  which  achieve  this  minimization. 

Nuttall  [10]  points  out  that  the  above  equation  is  minimized  by  setting  the 
coefficient  of  Gp  [contained  in  the  large  square  brackets  of  the  second  trace  term 

of  eq(4.55)]  equal  to  a  null  matrix  [15].  Premultiplying  this  coefficient  by  [Aplj] 

and  post  multiplying  by  [Qp-j],  we  obtain 


l^fblp-ilVp'-il[£Jp'-,l  +  Gp  [Vp'.,lH  (ftbblp-,lVp.il  £2p-i 

+  (Ap.,l[Up.,lH[ftn,|p.,  +  Ap.,[Up.,]H[ft„lp.,  [Up.,lGp  =  [0)  (4.56a) 

SO  that  after  rearranging  terms 

-  Gp  [Vp'.,]H  [ftbblp.,[Vp'-,l  £Jp‘.,  -  Ap-,|Up'.,]H[Aff)p.,  [Up'.,|Gp  = 

=  [^fbIp-i[Vp'-il[fip-il  +  [Ap'.,l|Up'.,]Hlftp,|p,,  ,  (4.56b) 

However,  from  eq(4.49) 

Gp  s  -  (Alp  (4.57) 

SO  that 

[A|pV  +  BH[&|«=cH  (4.58) 

A  H 

where  [A]p  is  recognized  as  an  estimate  since  limited  data  is  being  used  in  this 
equation  and 

A”  =  |Vp.,lH  [Abblp.,|Vp'.,)  Qp'.,  (4.59a) 

B«  =  Ap.,[Up.||H|ftff|p.,  [Up'.,|  (4.59b) 
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=  [ftfblp.|lV^l[£2p'-,l  +  [Ap‘-i)lUp-,]”|ftfblp-i  ■  (4.59c) 

In  reference  [12],  Nuttall  examined  several  forms  of  the  weighting  matrices 


■'^p-l 

and  Qp.j  and  showed  that  stable  conditions  are  achieved  if 

Ap.,  =  Up'.,slSf]p'., 

(4.60) 

and 

«p-i  =  Vp.,  s  (Z„|p'.,  . 

(4.61) 

Thus,  the  constant  matrices  A^,  B"  and  simplify  to 

A»  =  [Vp' ,|H  =  (Vp‘.,lH 

(4.62a) 

B”  =  (Afflp-i  [Up.,|  =  (Afflp.,  [Up.,]H 

(4.62b) 

C"  =  2[ftft,lp., 

(4.62c) 

where  [^bblp-f  [Up-il  are  Hermitian.  We  note  that  and  each 

contain  the  product  of  two  matrices.  The  matrices  [U’pl,]  and  [Vplj]  are  the 
forward  and  backward  error  variances  for  a  filter  of  order  p-1,  respectively. 
They  are  obtained  from  the  Levinson-Wiggins-Robinson  recursion.  The  matrices 
[ftfflp.j  and  [^bb^p-p  however,  are  time-averaged  estimates  of  these  variances 
using  limited  data  [see  eqs(4.47)].  Using  eqs(4.62),  eq(4.58)  is  now  written  as 

[^)“lVp'.,]H  [ftbbl"-!  +  [Afflp-,  [Up-i]"[A]p  =  (4.63) 

In  the  single  channel  case,  each  element  is  a  scalar  where 

Up.i  =  Vp.,  S  Op.,  (4.64) 


and  Ap  (p)  =  [Hflp  =  -  [A]^(Vp.|)H  becomes 

a;(p)  =  = 

^p-i 

From  eqs(4.47),  the  terms  [ftbblp-i » I^fflp-i  t^fblp-i 
expressed  as 


(4.65) 

for  the  scalar  case  can  be 
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(4.66a) 


*  1  N 

^  k=p+l 

A  1  N 

t^fblp-1  “  jsj.p  X  ^p-i(l^)Pp-i(^“l) 

^ k=p+l 

[^bblp-i  =  [ftbblp-i  =  I  IPp- ,  (k- 1  )l^- 

^ k=p+l 


(4.66b) 

(4.66c) 


Using  eqs(4.64)  through  (4.66),  the  scalar  version  of  eq(4.63)  can  be  written  as 


XiPp-i(k-l)|2 

- 

X  l£p-i 

a;(p)  =  2 

■  N 

£ep.i(k)Pp!i(k-l) 

_k=p+l 

k=p+l 

k=p+l 

(4.67) 

Solving  eq(4.67)  for  ap(p),  we  have 


N 

a^(P)  =  - •  (4.68) 

X[lep.,(k)|2  +  iPp.,(k-l)|2] 
k=p+l 

The  result  expressed  in  eq(4.68)  is  the  estimate  of  the  single  channel  Burg 

reflection  coefficient.  This  result  is  equivalent  to  that  obtained  by  minimizing 
eq(4.48b)  with  respect  to  Tp  [13].  Thus  the  results  presented  here  are  a 

generalization  of  the  single  channel  Burg  algorithm. 

4.4  SOME  CONSIDERATIONS  OF  BIASED  PARAMETER  ESTIMATORS 


In  this  section,  we  briefly  develop  a  relationship  that  will  be  useful  in 
evaluating  the  computed  estimation  results  in  chapter  5.0.  Specifically,  we 
develop  a  relationship  between  the  error  variance,  the  sample  variance  and  the 
bias  of  the  estimate. 

Consider  the  estimate  ^  of  the  parameter  a.  The  error  variance  of  ^  is 
expressed  as 
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V(,  =  E{[a-a][S-a]*}  (4.69a) 

=  E[iaP]  -  aE(a*)  -  a*E(a)  +  la|2  (4.69b) 

For  a  biased  estimator 

E(ft)  =  a  +  B(a)  (4.70a) 

E(a*)  =  a*  +  B*(a)  (4.70b) 

SO  that 

Vg  =  E[l^|2]  -  la|2  -  aB*(a)  -  a*B(a).  (4.71) 

Using  eqs(4.70),  the  variance  of  ^  is  expressed  as 

=  E{ [^  -  E(^  )][^  -  E(ll  )]* }  (4.72a) 

=  E[lll|2]  -  E(t)E(t*)  (4.72b) 

=  E[l^|2]  -  la|2  -  aB*(a)  -  a*B(a)  +lB(a)|2  (4.72c) 

Eqs(4.71)  and  (4.72c)  yield 

Vg  =  oj  +  IB(a)|2  (4.73) 


2 

The  quantity  can  be  estimated  by  the  computation  of  the  sample  variance 

of  In  chapter  5.0,  we  compute  the  quantities  in  eq(4.73)  for  verification  of  the 
results.  In  addition,  we  will  note  the  dependence  of  each  upon  the  correlation  of 
the  processes  for  various  estimators. 

4.5  SOME  ERROR  VARIANCE  CONSIDERATIONS  OF  PARAMETER 
ESTIMATORS 

As  noted  in  the  beginning  of  this  chapter,  expressions  for  the  error 
variance  of  time-averaged  parameter  estimators  are  difficult  to  determine.  This 
is  due  to  the  fact  that  the  estimators  generally  involve  ratios  of  correlated  random 
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estimates.  For  example,  from  the  Yule-Walker  equation,  the  estimate  ^i(l)  for 
an  AR(1)  process  is  expressed  as 


(4.74) 


while  for  an  AR(2)  process,  we  have  the  estimators 


and 


4  ...  _  ^(1)^(1)  -M^(0) 
^  ^(0)^(0)  -  ft(l)ft(-l) 


(4.75a) 


(4.75b) 


The  correlation  among  the  estimates  ift^(l)  is  due  to  the  use  of  a  single  set  of 
observation  data  in  the  estimation  procedure.  The  above  equations  reveal  a 
similar  form  to  that  noted  in  section  3.5  for  the  normalized  correlation  function 
estimators;  ie.,  the  deviation  of  the  numerator  terms  about  some  mean  level  will 
be  'tracked'  similarly  by  the  denominator  terms.  If  the  numerator  terms  increase 
(or  decrease),  the  denominator  terms  will  also  increase  (or  decrease).  As  a 
result,  a  type  of  'normalization'  results  so  that  the  estimator  will  be  closer  to  its 
respective  mean  and  its  associated  variance  is  reduced  as  compared  to  the 
corresponding  performance  of  the  individual  estimates  ft.(l). 
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5.0  RESULTS 


5.1  COMPUTED  VARIANCE  OF  THE  AUTOCORRELATION  FUNCTION 
ESTIMATOR 

In  this  section,  we  validate  eq(3.30)  which  is  a  special  case  of  eq(3.27) 
using  a  real,  exponentially  shaped  autocorrelation  function.  The  real  AR(1) 
process  also  has  an  exponential  autocorrelation  function  and  we  use  it  here  in  the 
synthesis  of  the  observation  processes.  Its  autocorrelation  function  is  expressed 
as 


RAR(k)  =  Rar(0)  [-ad)]"^' 

(5.1a) 

=  aiR  [-a(l)]"''. 

(5.1b) 

where 

1 

/.-w\  2 

=  1  -  a^(  1 )  ~ 

(5.2) 

2  2 

and  Ou  ,  a(l)  and  are  the  white  noise  driving  variance. 

the  AR(1)  parameter. 

and  the  variance  of  the  AR(1)  process,  respectively.  In  the  special  case  of  the 

AR(1)  process  used  in  this  example,  the  constants  in  eq(2.30)  were  chosen  such 
2  2 

that  A.jj=-a(l),  Oii=CFAR  and  9ii(l)=0  for  all  I.  Table  1  contains  the  parameters  used 

2 

in  the  process  synthesis  procedure  described  in  section  2.4.  The  variance  was 
held  fixed  at  4  while  X■^^  was  varied  using  0.1,  0.7  and  0.99. 


Fig. 

Nj 

Nr 

2 

4.0 

0.1 

100 

10,000 

3 

4.0 

0.7 

100 

10,000 

4 

4.0 

0.99 

100 

10,000 

Table  1  Parameters  used  in  the  synthesis  of  the  processes  analyzed  in  Figs.2 
through  4. 

Figs.  2  through  4  show  the  biased  time-averaged  autocorrelation  functions 
computed  with  eq(3.1)  and  their  associated  variances  for  the  processes  described 
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in  the  Table  1.  For  each  figure,  plot  (a)  shows  six  realizations  of  the  biased, 
time-averaged  autocorrelation  estimator  ^Tb(l)  plotted  over  64  lag  values  using 

N^=100  time  samples.  The  corresponding  ensemble  averaged  estimator  ^eCI)  is 
shown  in  plot  b  using  N(^=10,000.  The  sample  variance  Var[ltjj'p^(l,Nj):Nj^]  of  the 

biased  time-averaged  autocorrelation  function  computed  using  eq(3.28)  with 
Nr=10,000  is  displayed  in  c.  The  corresponding  analytic  calculation  of  Vb..(I,N-j-) 

using  eq(3.30)  is  shown  in  plot  d.  The  decrease  in  this  quantity  as  a  function  of  I 
is  a  result  of  the  weighting  associated  with  the  biased  estimator,  A  similar 
behavior  is  noted  in  [4]  for  continuous-time  processes.  The  computed  results 
(plots  c)  are  in  excellent  agreement  with  the  analytic  expressions  shown  in  plots  d 

(note  the  scale  change  between  plots  c  and  d  in  Fig.  2). 

In  Fig.  5,  the  maximum  value  of  Vr>  (I,N-t>)  which  occurs  at  1=0  is  plotted 

ii 

(solid  curve)  as  a  function  of  A,jj  for  Ny=100  and  Nj=1000  using  the  analytic 

expression  of  eq(3.30)  for  the  AR(1)  process.  The  corresponding  sample 
variances  of  the  time-averaged  autocorrelation  function  estimates  computed  at  lag 
zero  using  eq(3.28)  with  i=j  for  the  synthesized  data  processes  are  also  plotted  (•) 
on  this  curve.  These  values  were  computed  using  Nj^  realizations  of  the 
functions.  For  Nj  =1(X),  Nr=10,(X)0  was  used  while  for  Nj  =1000,  the  number 
of  realizations  was  reduced  to  Njj=l,000.  Reducing  Nj^  decreases  the  confidence 

level  associated  with  the  computation  of  the  sample  variance  Var[^ijj^(l,Ny):NR]. 

However,  the  error  bars  representing  one  standard  deviation  from  the  mean  are 
less  than  the  size  of  the  printed  (•), 

The  results  shown  in  Fig.5  indicate  the  significant  increase  in  the  variance 
Vr  (I,Nx)  as  a  function  of  X:,.  In  addition,  it  provides  a  measure  of  the  required 

“ii  * 

observation  window  size  necessary  to  achieve  a  specific  level  of  the  variance. 
For  example,  the  curve  for  N.j.=1000  has  a  very  distinct  knee  for  ^ji=0.9.  For 

processes  with  a  temporal  correlation  above  this  value,  a  larger  observation  data 
window  would  be  required  to  reduce  Vg.XI.Nj)  to  values  less  than  0.15.  We  also 

note  that  for  Nj=100,  large  variances  can  be  obtained  even  for  moderately  low 
values  of  In  [2],  similar  plots  are  shown  for  higher  values  of  Nj. 
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Figure  2  Time-averaged  autocorrelation  function  and  its  variance  for  an  AR(1) 
process;  A.=0.1,  cl=4  a.)  biased  (6  trials)  using  Nt=100  b.)  ensemble 

averaged  ^e(I)  using  10,000  realizations  c.)  computed  sample  variance  of  the 
biased  d.)  analytical  variance  of  the  biased 
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Figure  3  Time-averaged  autocorrelation  function  and  its  variance  for  an  AR(1) 
process;  >.=0.7,  aj=4  a.)  biased  (6  trials)  using  N.p=100  b.)  ensemble 

averaged  ^e(I)  using  10,000  realizations  c.)  computed  sample  variance  of  the 
biased  d.)  analytical  variance  of  the  biased 


Figure  4  Time-averaged  autocorrelation  function  and  its  variance  for  an  AR(1) 
process;  X,=0,99,  ^5^=4  a.)  biased  (4  trials)  using  Nj=100  b.)  ensemble 

averaged  ^e(I)  using  10,000  realizations  c.)  computed  sample  variance  of  the 
biased  d.)  analytical  variance  of  the  biased 


43 


Figure  5  Maximum  variance  of  the  time-averaged  autocorrelation  function 

2 

versus  for  ctjj=4;  analytical  (-)  and  computed  (•). 
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5.2  Computed  Variance  of  the  Cross-Correlation  Function 

Table  2  contains  the  parameters  used  in  the  process  synthesis  procedure 
[1,2]  for  a  two  channel  AR(1)  process  with  real  correlation  functions.  For  the 

cases  considered  here,  the  variances  of  the  two  processes,  denoted  as  On  and  O22’ 
were  both  held  fixed  at  4  while  and  had  values  of  0.1,  0.5  and  0.99.  The 
cross-correlation  coefficient  Ipj^l  had  values  of  0.99,  0.5  and  0.0. 


Fig. 

4 

^11 

^2 

^12 

'P.2' 

Nj 

Nr 

'12 

6 

4 

4 

0.1 

0.1 

0.1 

0.99 

100 

1000 

0 

7 

4 

4 

0.1 

0.1 

0.1 

0.50 

100 

1000 

0 

8 

4 

4 

0.1 

0.1 

0.1 

0.00 

100 

1000 

0 

9 

4 

4 

0.5 

0.5 

0.5 

0.50 

100 

1000 

0 

10 

4 

4 

0.9 

0.9 

0.9 

0.50 

100 

1000 

0 

Table  2  Parameters  used  in  the  synthesis  of  the  processes  analyzed  in  Figs.2 
through  4. 

Figs  6  through  10  show  the  variances  of  the  time-averaged  cross¬ 
correlation  function  estimates  based  on  the  computed  values  of  eq(3.28)  and  the 
analytic  expression  of  eq(3.31).  The  computed  results  are  in  excellent  agreement 
with  the  analytic  expressions.  In  Figs,6c,  7c  and  8c,  we  observe  that  the 
computed  sample  variance  is  not  affected  by  changes  in  Ipj2l.  This  obser\'ation  is 

consistent  with  the  point  made  in  section  3.2  that,  in  the  special  case  where 
eqs.(3.26)  hold,  Fjj(l,k)=0  and  the  variance  of  the  cross-correlation  function 

estimator  is  independent  of  the  cross-correlation  parameter  Ip  12!. 

In  Fig, 11,  the  maximum  value  of  Vg  (l,Nj),  which  occurs  at  I=lj2=0,  is 

plotted  (solid  curves)  as  a  function  of  X=A.,j=X,22  for  Nj=100  and  Nj=1000  using 

the  analytic  expression  of  eq(3.31)  and  a]j=G^2=^-  corresponding  maximum 
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:ross-coiTelation  function  and  its  variance  for  lo. 


•  c 
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Figure  10  Time-averaged  cross-correlation  function  and  its  variance  for  Ipj2l=0.50,  X,22=  ^i2=0,9, 

af, =022=4,  l|2=0  a.)  biased  (6  trials)  using  Np=100  b.)  ensemble  averaged  ^^(l)  using  1,000 

realizations  c.)  sample  variance  of  the  biased  d.)  analytical  variance  of  the  biased 


Figure  11  Variance  of  the  time-averaged  cross-correlation  function  at  l=l|2=0 

versus  X,=Xjj=X^2  ^22=4;  analytical(-)  and  computed  (•). 

values  of  the  sample  variances  of  the  time-averaged  cross-correlation  function 

estimates  computed  using  the  synthesized  data  processes  are  also  plotted  (•)  on 
this  curve.  These  values  were  computed  using  Nj^=  1000  realizations  of  the 

estimates.  In  several  of  the  figures,  we  note  the  increase  in  the  variance  of 

over  the  results  presented  in  the  previous  section.  This  is  simply  due  to  the 
decrease  in  from  10,000  to  1,000  thus  increasing  the  uncertainty  of  the 

statistical  calculations. 
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5.3  COMPUTED  ERROR  VARIANCE  OF  SINGLE  CHANNEL  AUTO¬ 
REGRESSIVE  MODEL  PARAMETERS 

In  this  section,  we  consider  complex  autoregressive  (AR)  time  series 
models  with  known  order  and  show  computed  simulation  results  for  the  error 
variance,  the  bias  and  the  variance  of  the  parameter  estimates  using  the  three 
estimators  described  in  chapter  4.0.  Both  Gaussian  and  K-distributed  processes 
are  considered. 

5.3.1  Performance  with  Gaussian  Processes 

Complex  Gaussian  autoregressive  (AR)  time  series  models  are  considered 
in  this  subsection.  Tables  3a,  b  and  c,  show  the  estimated  means,  error  variances 
and  sample  variances,  respectively,  of  the  a(l)  coefficient  for  Nj  values  ranging 

from  10  to  500.  These  values  were  obtained  for  each  estimator  using  a  complex 
single  channel  AR(1)  process  'without  noise'.  The  processes  were  synthesized 
with  the  real  exponential  correlation  function  of  eq(5.1b).  In  Tables  4a  through 
4f,  we  focus  on  N-j.=10  and  100  using  more  X  values.  The  value  of  the  bias  B  (the 

difference  between  the  estimates  ^(1)  and  the  true  value)  for  each  estimator  is 
also  shown.  The  method  described  in  section  2.4  was  used  to  synthesize  each 
process.  In  the  special  case  of  the  AR(1)  process,  the  one-lag  temporal 
correlation  parameter  A,  is  the  negative  of  the  AR  coefficient.  Although  the  true 
parameters  of  the  autoregressive  process  are  real,  the  processes  themselves  are 
complex.  Thus,  an  imaginary  component  of  the  estimate  will,  in  general,  be 
obtained.  This  component,  however,  becomes  negligible  in  the  limit  as  the 
estimate  approaches  the  true  value.  Only  the  means  of  the  real  part  of  the 
estimated  coefficients  are  shown  here  since  the  complex  parts  were  small  for 
these  cases.  Computed  values  of  the  error  variance  of  ^(1)  and  Ay  are  obtained 
using  N(j=  10,000  realizations  where  S(l)  and  6,,  are  the  estimates  of  a(l)  and  0^, 
respectively.  Fig.  12a  shows  the  error  variances  of  S(l)  versus  A  using  Nj=10  for 
the  three  estimators.  Corresponding  log  scale  plots  for  Ny=10  and  100  are  shown 

in  Figs.  12b  and  12c,  respectively.  Each  of  these  plots  are  compared  to  the  exact 
unbiased  Cramer-Rao  bound  [6].  We  note  that  the  error  variances  of  the 
estimators  lie  below  that  of  the  Cramer-Rao  bound  for  most  values  of  A.  This 
observation  is  not  unexpected  since  these  estimators  are  not  unbiased. 
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0.9999 


a 


-  0.1 


-  0.1 


-  0.1 


-  0.1 


-  0.1 


-0.5 


-0.5 


-  0.5 


-  0.5 


-  0.5 


-  0.9 


-  0.9 


-0.9 


-0.9 


-0.9 


-  0.99 


-  0.99 


-0.99 


-  0.99 


-  0.99 


0.9999 


0.9999 


0.9999 


0.9999 


0.9999 


Nj 

mean  Re  5(1) 
YWUBC 

mean  Re  5(1) 
YWBC 

10 

-  8.9437x10-2 

-  8.0494x10-2 

50 

-  9.9174x10-2 

-  9.7190x10-2 

100 

-  9.9981x10-2 

-9.8981x10-2 

200 

-  1.0006x10-1 

-  9.9560x10-2 

500 

-  1.0029x10-1 

-  1.0009x10-1 

10 

-  4.5493x10-1 

-4.0944x10-1 

50 

-  4.9035x10-1 

-  4.8055x10-1 

100 

-  4.9533x10-1 

-  4.9038x10-1 

200 

-  4.9724x10-1 

-  4.9475x10-1 

500 

-  4.9872x10-1 

-  4.9773x10-1 

10 

-  8.3169x10-1 

-  7.4852x10-1 

50 

-  8.8284x10-1 

-8.6518x10-1 

100 

-  8.9088x10-1 

-  8.8197x10-1 

200 

-  8.9557x10-1 

-  8.9109x10-1 

500 

-  8.9802x10-1 

-  8.9622x10-1 

10 

-  9.6398x10-1 

-  8.6758x10-1 

50 

-  9.7772x10-1 

-  9.5817x10-1 

100 

-  9.8209x10-1 

-  9.7226x10-1 

200 

-  9.8573x10-1 

-  9.8080x10-1 

500 

-  9.8801x10-1 

-  9.8603x10-1 

10 

-  9.9896x10-1 

-  8.9906x10-1 

50 

-  9.9929x10-1 

-9.7931x10-1 

100 

-  9.9931x10-1 

-  9.8932x10-1 

200 

-  9.9945x10-1 

-  9.9446x10-1 

500 

-  9.9955x10-1 

-  9.9755x10-1 

-  8.9294x10-2 


-  9.9191x10-2 


-  1.0031x10-1 


-  1.0006x10-1 


-  1.0029x10-1 


-  4.5390x10-1 


-  4.9029x10-1 


-  4.9534x10-1 


-4.9722x10-1 


-  4.9873x10-1 


-  8.3815x10-1 


-  8.8369x10-1 


-8.9110x10-1 


-  8.9562x10-1 


-  8.9802x10-1 


-  9.6822x10-1 


-  9.7979x10-1 


-  9.8332x10-1 


-  9.8623x10-1 


-9.8813x10-1 


-  9.9917x10-1 


-  9.9947x10-1 


-  9.9949x10-1 


-  9.9959x10-1 


-  9.9967x10-1 


Table  3a  Tabulated  values  of  the  mean  for  the  real  part  of  S(l)  with  X  and  as 
parameters  and  computed  using  Nj^=l  0,000  realizations. 
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X 

Nj 

error  variance 

ofS(l) 

YWUBC 

error  variance 
OfS(l) 
YWBC 

error  variance 
ofS(l) 

Burg  Algorith. 

unbiased  exact 

Cramer-Rao 

bound 

0.1 

10 

9.9742x10-2 

8.1081x10-2 

9.9524x10-2 

1.1088x10-1 

50 

1.94440x10-2 

1.8678x10-2 

1.9449x10-2 

100 

9.9574x10-3 

9.7603x10-3 

9.8274x10-3 

1.0099x10-2 

200 

4.9522x10-3 

4.9030x10-3 

4.9518x10-3 

500 

1.9872x10-3 

1.9791x10-3 

1.9871x10-3 

0.5 

10 

8.4426x10-2 

7.4780x10-2 

8.2327x10-2 

1.04167x10-1 

50 

1.5339x10-2 

1.5021x10-2 

1.5196x10-2 

100 

7.7212x10-3 

7.6388x10-3 

7.6932x10-3 

1.0034x10-2 

200 

3.7136x10-3 

3.6965x10-3 

3.7058x10-3 

500 

1.5107x10-3 

1.5082x10-3 

1.5088x10-3 

0.9 

10 

4.0044x10-2 

5.1605x10-2 

3.5420x10-2 

5.9975x10-2 

50 

5.1584x10-3 

5.8837x10-3 

4.7416x10-3 

100 

2.2807x10-3 

2.4789x10-3 

2.1737x10-3 

9.3074x10-3 

200 

1.0532x10-3 

1.1026x10-3 

1.0269x10-3 

500 

3.9503x10-4 

4.0380x10-4 

3.8968x10-4 

0.99 

10 

8.2635x10-3 

2.1132x10-2 

6.6868x10-3 

1.024x10-2 

50 

1.1963x10-3 

2.0175x10-3 

9.1984x10-4 

100 

4.8610x10-4 

7.2965x10-4 

3.8344x10-4 

5.0089x10-3 

200 

1.8450x10-4 

2.4927x10-4 

1.5359x10-4 

500 

5.7238x10-5 

6.8802x10-5 

5.1501x10-5 

0.9999 

10 

1.7938x10-4 

1.0314x10-2 

1.3863x10-4 

1.11x10-4 

50 

2.4836x10-5 

4.4756x10-4 

1.7980x10-5 

100 

1.5377x10-5 

1.2665x10-4 

1.0509x10-5 

1.0x10-4 

200 

5.3456x10-6 

3.4784x10-5 

3.8399x10-6 

500 

2.2383x10-6 

7.6558x10-6 

1.6141x10-6 

Table  3b  Tabulated  values  of  the  error  variances  of  S(l)  with  X  and  Nj  as 
parameters  and  computed  using  Nj^=10,000  realizations. 
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0.9999 


Nj 

1 

sample  variance 
ofS(l) 
YWUBC 

sample  variance 
ofS(l) 
YWBC 

sample  variance 
ofS(l) 

Burg  Algorith. 

10 

9.9627x10-2 

8.0698x10-2 

9.9406x10-2 

50 

1.9437x10-2 

1.8667x10-2 

1.9446x10-2 

100 

9.9569x10-3 

9.7588x10-3 

9.8249x10-3 

200 

4.9521x10-3 

4.9027x10-3 

4.9516x10-3 

500 

1.9871x10-3 

1.9791x10-3 

1.9870x10-3 

10 

8.2194x10-2 

6.6577x10-2 

8.0201x10-2 

50 

1.5246x10-2 

1.4642x10-2 

1.5102x10-2 

100 

7.6993x10-3 

7.5461x10-3 

7.6714x10-3 

200 

3.7043x10-3 

3.6674x10-3 

3.6965x10-3 

500 

1.5090x10-3 

1.5029x10-3 

1.5070x10-3 

10 

3.5373x10-2 

2.8652x10-2 

3.1590x10-2 

50 

4.8637x10-3 

4.6711x10-3 

4.4756x10-3 

100 

2.1974x10-3 

2.1537x10-3 

2.0945x10-3 

200 

1.0336x10-3 

1.0233x10-3 

1.0077x10-3 

500 

3.9107x10-4 

3.8951x10-4 

3.8572x10-4 

10 

7.5861x10-3 

6.1448x10-3 

6.2120x10-3 

50 

1.0455x10-3 

1.0041x10-3 

8.1546x10-4 

100 

4.2336x10-4 

4.1493x10-4 

3.3881x10-4 

200 

1.6623x10-4 

1.6457x10-4 

1.3940x10-4 

500 

5.3270x10-5 

5.3057x10-5 

4.7999x10-5 

10 

1.7848x10-4 

1.4457x10-4 

1.3803x10-4 

50 

2.4469x10-5 

2.3500x10-5 

1.7759x10-5 

100 

1.5035x10-5 

1.4735x10-5 

1.0299x10-5 

200 

5.1416x10-6 

5.0902x10-6 

3.7203x10-6 

500 

2.1105x10-6 

2.1021x10-6 

1.5371x10-6 

Table  3c  Tabulated  values  of  the  sample  variances  of  S(l)  with  X  and  as 
parameters  and  computed  using  N,^= 10,000  realizations. 
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Nj 

true 

a(l) 

mean 

ReS(l) 

YWBC 

error  variance 
ofS(l) 
YWBC 

10 

-0.1 

-8.0494x10-2 

8.1081x10-2 

100 

-0.1 

-9.8981x10-2 

9.7603x10-3 

10 

-0.2 

-1.6318x10-1 

8.0896x10-2 

100 

-0.2 

-1.9633x10-1 

9.6344x10-3 

10 

-0.3 

-2.4547x10-1 

7.8376x10-2 

100 

-0.3 

-2.9440x10-1 

9.0049x10-3 

10 

-0.4 

-3.2317x10-1 

7.7914x10-2 

100 

-0.4 

-3.9142x10-1 

8.3293x10-3 

10 

-0.5 

-4.0944x10-1 

7.4780x10-2 

100 

-0.5 

-4.9533x10-1 

7.6388x10-3 

10 

-0.6 

-4.8869x10-1 

7.3030x10-2 

100 

-0.6 

-5.8840x10-1 

6.5236x10-3 

10 

-0.7 

-5.7702x10-1 

6.6399x10-2 

100 

1 

o 

-6.8598x10-1 

5.3763x10-3 

10 

-0.8 

-6.6043x10-1 

6.1667x10-2 

100 

-0.8 

-7.8483x10-1 

3.9796x10-3 

10 

-0.85 

-7.0213x10-1 

5.6941x10-2 

100 

-0.85 

-8.3350x10-1 

3.2694x10-3 

10 

-0.90 

-7.4852x10-1 

5.1605x10-2 

100 

-0.90 

-8.8197x10-1 

2.4789x10-3 

10 

-0.95 

-8.0363x10-1 

4.0403x10-2 

100 

-0.95 

-9.3093x10-1 

1.6403x10-3 

10 

-0.99 

-8.6758x10-1 

2.1132x10-2 

100 

-0.99 

-9.7226x10-1 

7.2965x10-^ 

10 

-0.9999 

-8.9906x10-1 

'  1.0314x10-2 

100 

-0.9999 

-9.8932x10-1 

1.2665x10-^ 

sample  variance 
ofS(l) 
YWBC 


8.0698x10-2 


9.7588x10'^ 


7.9548x10-2 


9.6219x10-3 


7.5409x10-2 


8.9737x10-3 


7.2017x10-2 


8.2560x10-3 


6.6577x10-2 


7.5461x10-2 


6.0644x10-2 


6.3896x10-3 


5.1278x10-2 


5.1804x10-3 


4.2189x10-2 


3.7496x10-3 


3.5080x10-2 


2.9930x10-3 


2.8652x10-2 


2.1537x10-3 


1.8980x10-2 


1.2766x10-3 


6.1448x10-3 


1.4457x10-"^ 


1.4735x10-3 


1.9506x10-2 


1.0190x10-3 


3.6820x10-2 


3.6700x10-2 


5.4530x10-2 


5.6000x10-2 


7.6830x10-2 


8.5800x10-3 


9.0560x10-2 


4.6700x10-3 


1.1131x10-* 


1.1600x10-2 


1.2298x10-1 


1.4020x10-2 


1.3957x10-1 


1.5170x10-2 


1.4787x10-1 


1.6500x10-2 


1.5148x10-1 


1.8030x10-2 


1.4637x10-1 


1.9070x10-2 


1.2242x10-1 


1.7740x10-2 


1.0084x10-1 


1.0580x10-2 


Table  4a  Tabulated  values  of  the  statistics  of  S(l)  for  the  YWBC  with  X.  and  Nj  as 
parameters  and  computed  using  Nr=  10,000  realizations. 
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X 

Nj 

true 

mean 

YWBC 

error  variance 
of^ 

YWBC 

sample  variance 
of  6^ 

YWBC 

bias 

YWBC 

0.1 

10 

3.96 

3.6359 

1.5455 

1.4405 

3.2410x10'* 

100 

3.96 

3.9215 

1.6015x10* 

1.5867x10-* 

3.8500x10-2 

0.2 

10 

3.84 

3.5614 

1.4756 

1.3981 

2.7860x10'* 

100 

3.84 

3.8092 

1.4627x10-* 

1.4534x10-* 

3.0800x10-2 

0.3 

10 

3.64 

3.4152 

1.3107 

1.2603 

2.2480x10-* 

100 

3.64 

3.6099 

1.3955x10  * 

1.3865x10'* 

3.0100x10-* 

0.4 

10 

3.36 

3.1971 

1.1324 

1.1060 

1.6290x10-* 

100 

3.36 

3.3389 

1.1306x10-* 

1.1262x10-* 

2.1100x10-2 

0.5 

10 

3.00 

2.9463 

9.7856x10-* 

9.7567x10-* 

5.3700x10-2 

100 

3.00 

2.9920 

8.8970x10-2 

8.8906x10-2 

8.0000x10-3 

0.6 

10 

2.56 

2.6140 

7.7840x10-* 

7.7556x10-* 

-5.4000x10-2 

100 

2.56 

2.5659 

6.7719x10-2 

6.7691x10-2 

-5.9000x10-3 

0.7 

10 

2.04 

2.2396 

6.5311x10-* 

6.1332x10-* 

-1.9960x10'* 

100 

2.04 

2.0578 

4.4638x10-2 

4.4324x10-2 

-1.7800x10-2 

0.8 

10 

1.44 

1.797 

5.9205x10-* 

4.6463x10-* 

-3.5700x10'* 

100 

1.44 

1.4776 

2.4925x10-2 

2.3512x10-2 

-3.7600x10-2 

0.85 

10 

1.11 

1.5772 

6.3943x10-* 

4.2122x10'* 

-4.6720x10'* 

100 

1.11 

1.1602 

1.8008x10-2 

1.5490x10-2 

-5.0200x10-2 

0.9 

10 

0.76 

1.3101 

7.0149x10-* 

3.9884x10'* 

-5.5010x10* 

100 

0.76 

8.1552x10-^ 

8.7135x10-^ 

1.1796x10-2 

-5.5520x10-2 

0.95 

10 

0.39 

1.0394 

8.3888x10-* 

4.1722x10'* 

-6.4940x10-* 

100 

0.39 

4.5763x10-* 

9.2025x10-3 

4.6284x10-3 

-6.7630x10-2 

0.99 

10 

0.0796 

8.2175x10-* 

1.0869 

5.3602x10  * 

-7.4215x10-* 

100 

0.0796 

1.5669x10* 

9.6108x10-3 

3.6671x10-3 

-7.7090x10-2 

0.9999 

10 

8.0013x-^ 

7.5213x10-* 

1.1240 

5.5949x10* 

-7.5133x10'* 

100 

8.0013x-^ 

8.0369x10-2 

1.2688x10-2 

6.3565x10-3 

-7.9569x10-2 

Table  4b  Tabulated  values  of  the  statistics  of  6^  for  the  YWBC  with  X  and  Nj  as 
parameters  and  computed  using  Nr= 10,000  realizations. 


Nt 

H 

mean 

Re^(l) 

YWUBC 

error  variance 

ofa(l) 

YWUBC 

sample  variance 
ofS(l) 

YWUBC 

bias 

B[S(1)| 

YWUBC 

0.1 

10 

-0.1 

-8.9437x10-2 

9.9742x10-2 

9.9627x10-2 

1.0563x10-2 

100 

-0.1 

-9.9981x10-2 

9.9574x10-2 

9.9569x10-2 

1.9000x10-5 

0.2 

10 

-0.2 

-1.8161x10-* 

9.9547x10-2 

9.9218x10-2 

1.8390x10-2 

100 

-0.2 

-1.9781x10-* 

9.6865x10-2 

9.6827x10-2 

2.1900x10-2 

0.3 

10 

-0.3 

-2.7341x10-* 

9.5300x10-2 

9.4582x10-2 

2.6590x10-2 

100 

-0.3 

-2.9687x10-* 

9.0241x10-2 

9.0143x10-2 

3.1300x10-2 

0.4 

10 

-0.4 

-3.6434x10-* 

9.0698x10-2 

8.9435x10-2 

3.5660x10-2 

100 

-0.4 

-3.9608x10  * 

8.6024x10-2 

8.5877x10-2 

3.9200x10-2 

0.5 

10 

-0.5 

-4.5493x10-* 

8.4426x10-2 

8.2194x10-2 

4.5070x10-2 

100 

-0.5 

-4.9533x10-* 

7.7212x10-2 

7.6993x10-2 

4.6700x10-2 

0.6 

10 

-0.6 

-5.4432x10-* 

7.7740x10-2 

7.4647x10-2 

5.5680x10-2 

100 

-0.6 

-5.9340x10-* 

6.5662x10-2 

6.5231x10-2 

6.6000x10-2 

0.7 

10 

-0.7 

-6.3791x10  * 

6.7425x10-2 

6.3568x10-2 

6.2090x10-2 

100 

-0.7 

-6.9267x10-* 

5.2886x10-2 

5.2353x10-2 

7.3300x10-2 

0.8 

10 

-0.8 

-7.3311x10-* 

5.5550x10-2 

5.1078x10-2 

6.6890x10-2 

100 

-0.8 

-7.9232x10-* 

3.9145x10-2 

3.8551x10-2 

7.6800x10-2 

0.85 

10 

-0.85 

-7.8411x10-* 

4.7431x10-2 

4.3093x10-2 

6.5890x10-2 

100 

-0.85 

-8.4202x10-* 

3.1096x10-2 

3.0453x10-2 

7.9800x10-2 

0.90 

10 

-0.9 

-8.3169x10* 

4.0044x10-2 

3.5373x10-2 

6.8310x10-2 

100 

-0.9 

-8.9088x10-* 

2.2807x10-2 

2.1974x10-2 

9.1200x10-2 

0.95 

10 

-0.95 

-8.9609x10* 

2.5219x10-2 

2.2315x10-2 

5.391x10-2 

100 

-0.95 

-9.4101x10* 

1.3332x10-2 

1.2522x10-2 

8.9900x10-2 

0.99 

10 

-0.99 

-9.6398x10* 

8.2635x10-2 

7.5861x10-2 

2.6020x10-2 

100 

-0.99 

-9.8209x10-* 

4.8610x10-^ 

4.2336x10-4 

7.9100x10-2 

0.9999 

10 

-0.9999 

-9.9896x10* 

1.7938x10-4 

1.7848x10-4 

9.4000x10-4 

100 

-0.9999 

-9.9931x10* 

1.5377x10-2 

1.5035x10-^ 

5.9000x10-4 

Table  4c  Tabulated  values  of  the  statistics  of  ^(1)  for  the  YWUBC  with  X  and 
as  parameters  and  computed  using  Njj=  10,000  realizations. 
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X 

true 

mean  Re  6^ 

YWUBC 

error  variance 
OfftJ 

YWUBC 

sample  variance 
OfftJ 

YWUBC 

bias 

B[f^J 

YWUBC 

0.1 

10 

3.96 

3.5536 

1.5976 

1.4325 

4.0640x10*’ 

100 

3.96 

3.9199 

1.6020x10-’ 

1.5860x10’ 

4.0100x10-2 

0.2 

10 

3.84 

3.4482 

1.4830 

1.3297 

3.9180x10’ 

100 

3.84 

3.8051 

1.4842x10-’ 

1.4722x10-’ 

3.4900x10*2 

0.3 

10 

3.64 

3.2627 

1.3887 

1.2465 

3.7730x10*’ 

100 

3.64 

3.6099 

1.2943x10-’ 

1.2853x10-’ 

3.0100x10*2 

0.4 

10 

3.36 

3.0185 

1.1599 

1.0433 

3.4150x10*’ 

100 

3.36 

3.3219 

1.1348x10-’ 

1.1205x10-’ 

3.8100x10-2 

0.5 

10 

3.00 

2.6973 

9.6695x10’ 

8.7529x10-’ 

3.0270x10  ’ 

100 

3.00 

2.9714 

8.8579x10-2 

8.7761x10-2 

2.8600x10*2 

0.6 

10 

2.56 

2.2974 

7.2121x10-’ 

6.5231x10-’ 

2.6260x10*’ 

100 

2.56 

2.5329 

6.5913x10-2 

6.5186x10-2 

2.7100x10-2 

0.7 

10 

2.04 

1.8303 

5.2519x10-’ 

4.8125x10*’ 

2.0970x10"’ 

100 

2.04 

2.0199 

4.2119x10-2 

4.1721x10-2 

2.0100x10-2 

0.8 

10 

1.44 

1.2782 

3.4662x10-’ 

3.2047x10-’ 

1.6180x10*’ 

100 

1.44 

1.4227 

2.1857x10-’ 

2.1560x10*2 

1.7300x10*2 

0.85 

10 

1.11 

0.98478 

2.6197x10-’ 

2.4632x10*’ 

1.2522x10-’ 

100 

1.11 

1.0950 

1.3562x10-2 

1.3349x10*2 

1.5000x10*2 

0.9 

10 

0.76 

6.8370x10-’ 

1.8335x10-’ 

1.7752x10*’ 

7.6300x10-2 

100 

0.76 

7.5113x10-’ 

7.8289x10-3 

7.7503x10*3 

8.8700x10-3 

0.95 

10 

0.39 

3.5411x10-’ 

9.9623x10-2 

9.8345x10*2 

3.5890x10-2 

100 

0.39 

3.8545x10-’ 

3.9760x10-3 

3.9557x10*3 

4.5500x10-3 

0.99 

10 

0.0796 

7.3049x10-2 

1.8724x10-2 

1.8681x10*2 

6.5510x10-3 

100 

0.0796 

7.8683x10-2 

1.7837x10-3 

1.7828x10*3 

9.1700x10*4 

0.9999 

10 

8.001 3x10''^ 

6.6967x10-^ 

1.8188x10-4 

1.8186x10*4 

1.3046x10-4 

100 

8.001 3x10-'^ 

7.8370x10-4 

2.1419x10-3 

2.1419x10*3 

1.6430x10-3 

Table  4d  Tabulated  values  of  the  statistics  of  for  the  YWUBC  with  X  and  Nj  as 
parameters  and  computed  using  Nr= 10,000  realizations. 
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X 

Nj 

true 

a(l) 

mean 

ReS(l) 

BURG 

error  variance 
of&(l) 
BURG 

sample  variance 
ofS(l) 
BURG 

bias 

B[9(l)l 

BURG 

0.1 

10 

-0.1 

-8.9294x10-2 

9.9524x10-2 

9.9406x10*2 

1.0706x10*2 

100 

-0.1 

-1.0031x10-* 

9.8274x10-2 

9.8249x10*2 

-3.1000x10*4 

0.2 

10 

-0.2 

-1.8327x10-* 

9.9120x10-2 

9.8850x10*2 

1.6730x10*2 

100 

-0.2 

-1.9729x10-* 

9.6306x10-2 

9.6242x10*2 

2.7100x10*2 

0.3 

10 

-0.3 

-2.7314x10-* 

9.3171x10-2 

9.2458x10*2 

2.6860x10*2 

100 

-0.3 

-2.9661x10-* 

9.2899x10-2 

9.2793x10*2 

3.3900x10*2 

0.4 

10 

-0.4 

-3.6327x10-* 

8.8534x10-2 

8.7192x10*2 

3.6730x10*2 

100 

-0.4 

-3.9625x10-* 

8.4301x10*2 

8.4168x10*2 

3.7500x10*2 

0.5 

10 

-0.5 

-4.5390x10  * 

8.2327x10-2 

8.0201x10*2 

4.6100x10*2 

100 

-0.5 

-4.9534x10-* 

7.6932x10-2 

7.6714x10-2 

4.6600x10*2 

0.6 

10 

-0.6 

-5.4378x10-* 

7.5013x10-2 

7.1847x10*2 

5.6220x10-2 

100 

-0.6 

-5.9440x10-* 

6.4842x10-2 

6.4531x10-2 

5.6000x10*2 

0.7 

10 

-0.7 

-6.3707x10-* 

6.4977x10-2 

6.1022x10-2 

6.2930x10*2 

100 

-0.7 

-6.9288x10-* 

5.2134x10-2 

5.1632x10*2 

7.1200x10*2 

0.8 

10 

-0.8 

-7.3561x10  * 

5.1247x10-2 

4.7091x10*2 

6.4390x10*2 

100 

-0.8 

-7.9190x10-* 

3.8491x10-2 

3.7838x10*2 

8.1000x10*2 

0.85 

10 

-0.85 

-7.8489x10-* 

4.3703x10*2 

3.9469x10*2 

6.5110x10*2 

100 

-0.85 

-8.4187x10-* 

3.0015x10*2 

2.9356x10*2 

8.1300x10*2 

0.9 

10 

-0.9 

-8.3815x10-* 

3.5420x10*2 

3.1590x10*2 

6.1850x10*2 

100 

-0.9 

-8.9110x10-* 

2.1737x10*2 

2.0945x10*2 

8.9000x10*2 

0.95 

10 

-0.95 

-8.9749x10-* 

2.3537x10*2 

2.0781x10*2 

5.2510x10*2 

100 

-0.95 

-9.4120x10-* 

1.2599x10-2 

1.1822x10*2 

8.8000x10*2 

0.99 

10 

-0.99 

-9.6822x10-* 

6.6868x10*2 

6.2120x10-2 

2.1780x10*2 

100 

-0.99 

-9.8332x10* 

3.8344x10'^ 

3.3881x10*^ 

6.6800x10*2 

0.9999 

10 

-0.9999 

-9.9917x10* 

1.3863x10-^ 

1.3803x10*4 

7.3000x10*4 

100 

-0.9999 

-9.9949x10-* 

1.0509x10*5 

1.0299x10*5 

4.1000x10*4 

Table  4e  Tabulated  values  of  the  statistics  of  ^(1)  for  the  Burg  estimator  with  X 
and  Ny  as  parameters  and  computed  using  Nj^=  10,000  realizations. 
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X 

Nj 

true 

mean  Re  6^ 

BURG 

error  variance 
of^ 

BURG 

sample  variance 
of  6^ 

BURG 

bias 

BURG 

0.1 

10 

3.96 

3.5694 

1.6092 

1.4566 

3.9060x10  * 

100 

3.96 

3.9199 

1.6018x10'* 

1.5858x10'* 

4.0100x10-2 

0.2 

10 

3.84 

3.4634 

1.4662 

1.3245 

3.7660x10  * 

100 

3.84 

3.7938 

1.4626x10'* 

1.4414x10'* 

4.6200x10-2 

0.3 

10 

3.64 

3.2607 

1.3007 

1.1569 

3.7930x10* 

100 

3.64 

3.6014 

1.3610x10* 

1.3463x10'* 

3.8600x10-2 

0.4 

10 

3.36 

3.0320 

1.1482 

1.0407 

3.2800x10  * 

100 

3.36 

3.3264 

1.1364x10'* 

1.1252x10'* 

3.3600x10-2 

0.5 

10 

3.00 

2.7151 

9.1918x10* 

8.3800x10'* 

2.8490x10'* 

100 

3.00 

2.9715 

8.8333x10'2 

8.7523x10-2 

2.8500x10-2 

0.6 

10 

2.56 

2.3152 

6.5754x10'* 

5.9765x10'* 

2.4480x10'* 

100 

2.56 

2.5322 

6.4634x10-2 

6.3870x10-2 

2.7800x10-2 

0.7 

10 

2.04 

1.8583 

4.1647x10'* 

3.8349x10'* 

1.8170x10'* 

100 

2.04 

2.0216 

5.2134x10-3 

4.0778x10-2 

1.8400x10-2 

0.8 

10 

1.44 

1.3143 

2.1044x10'* 

1.9466x10-* 

1.2570x10'* 

100 

1.44 

1.4238 

2.0554x10-2 

2.0294x10-2 

1.6200x10-2 

0.85 

10 

1.11 

1.0187 

1.2332x10'* 

1.1500x10'* 

9.1300x10-2 

100 

1.11 

1.0992 

1.2353x10-2 

1.2238x10-2 

8.0000x10-2 

0.9 

10 

0.76 

6.9803x10'* 

5.9288x10-2 

5.5448x10-2 

6.1970x10-2 

100 

0.76 

7.5200x10'* 

5.7107x10-3 

5.6467x10-3 

8.0000x10-3 

0.95 

10 

0.39 

3.6265x10'* 

1.5734x10-2 

1.4988x10-2 

2.735x10-2 

100 

0.39 

3.8643x10'* 

1.5107x10-3 

1.4981x10-3 

3.5700x10-3 

0.99 

10 

0.0796 

7.4523x10'2 

6.5985x10'^ 

6.3407x1 0-4 

5.0770x10-3 

100 

0.0796 

7.9082x1 0'2 

6.4060x10-5 

6.3791x10-5 

5.1800x10-^ 

0.9999 

10 

S.OOOx''^ 

7.5 1  84x1  0'^ 

6.9146x10-8 

6.6810x10-8 

4.8290x10-5 

100 

S-OOnx'"* 

7.9691x10'^ 

6.3769x10-9 

6.3663x10-9 

3.2200x10-2 

Table  4f  Tabulated  values  of  the  statistics  of  6^  for  the  Burg  estimator  with  X  and 
Nj  as  parameters  and  computed  using  Nj^=  10,000  realizations. 
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ERROR  VARIANCE 
OF  A(l)  ESTIMATE 


Figure  12a  Error  variance  of  ^l(l)  for  the  time  averaged  parameter  estimators  of 
an  AR(1)  process  versus  the  one-lag  temporal  correlation  parameter 
X  using  Nj=10. 
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ERROR  VARIANCE 
OF  A(l)  ESTIMATE 


YWUBC 

YWBC 

BURG 

EXACT  UNBIASED  CR  BOUND 


Figure  12b  Error  variance  of  ll(l)  for  the  time  averaged  parameter  estimators  of 
an  AR(1)  process  versus  the  one-lag  temporal  correlation  parameter 
X,  using  Nj=10  plotted  on  a  log  scale. 


ERROR  VARIANCE 
OF  A(l)  ESTIMATE 


Figure  12c  Error  variance  of  ^(1)  for  the  time  averaged  parameter  estimators  of 
an  AR(1)  process  versus  the  one-lag  temporal  correlation  parameter 
X  using  N-p=100. 
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Unlike  the  correlation  function  estimators,  the  error  variance  of  the  AR 
parameter  estimators  decreases  with  increasing  temporal  correlation.  Between 
X,=0.9  and  0.9999,  the  decrease  in  error  variance  for  these  estimators  is  nearly 

three  orders  of  magnitude  as  shown  by  the  log  scale  plots  of  Figs.  12b  and  12c.  In 
Fig, 12a,  it  is  also  interesting  to  note  that  for  the  small  Nj  size  used  here  (ie., 
Nj=10),  the  error  variance  for  (YWBC)  is  the  lowest  of  the  three  for  processes 

with  low  temporal  correlation  (ie.,  X~  0.1).  However,  in  Fig. 13,  we  also 
observe  that  the  bias  B[^(l)]  of  the  (YWBC)  estimate  5(1)  computed  over  the 
10,000  realizations  is  larger  compared  to  the  other  two  estimators  for  all  values 
of  X  and  Nj.  The  explanation  for  this  apparently  contradictory  behavior  is  found 

in  Fig.  14  where  we  observe  that  the  variance  of  5(1)  using  (YWBC)  is  much 

smaller  than  the  other  estimators.  From  eq(4.73)  in  section  4.4,  we  note  that  a 

2 

very  small  variance  can  compensate  for  a  large  bias  thus  resulting  in  a  low 

error  variance.  At  a  value  of  ^  =  0.65,  however,  a  cross-over  occurs  in  the 

curves  shown  in  Fig.  12a  so  that  the  error  variances  of  the  Burg  and  (YWUBC) 

algorithms  are  lower  than  the  (YWBC)  at  higher  temporal  correlation  values. 

This  results  from  the  lower  biases  of  the  former  algorithms  as  well  as  the  fact 
2 

that  the  values  for  the  three  estimators  converge  at  high  X  (see  Fig.  14). 

In  Figs.  15a,  b  and  c,  we  show  the  convergence  of  the  mean  estimate  5(1)  to 
the  true  parameter  value  as  N-p  increases  for  values  of  A=0.1,  0.9,  and  0.9999, 
respectively.  For  small  time  window  sizes  (ie.,  small  Ny),  the  bias  of  the 

YWUBC  is  barely,  but  slightly,  smaller  than  that  of  the  Burg  estimator  (see 

Fig.  15a).  As  the  temporal  correlation  increases,  however,  the  performance  of  the 

Burg  estimator  is  superior  as  shown  in  Figs.  15b  and  c.  Finally,  although  the  bias 

of  the  YWBC  estimator  appears  to  be  increasing  dramatically  for  increasing 

temporal  correlation  (see  Fig.  15c)  as  compared  to  the  other  estimators,  the 

percentage  error  is  actually  decreasing  thus  showing  performance  improvement 

with  increasing  temporal  correlation.  In  Figs.  16,  17  and  18,  we  show  plots  for 

the  error  variance,  the  bias  and  the  estimator  variance  for  6^. 

Tables  5a,  b  and  c  list  the  means,  error  variances  and  estimator  variances 
2  2 

,  respectively,  of  5(1),  5(2)  and  Ay  for  single  channel  AR(2)  processes  with 
various  one-lag  temporal  correlation  parameters  X  using  the  order  2  Burg 
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X 

Figure  14a  Sample  variance  of  ^(1)  for  the  time  averaged  parameter  estimators 

of  an  AR(1)  process  versus  the  one-lag  temporal  correlation 
parameter  X  using  Nj=10. 
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Figure  15c  Mean  of  S(l)  versus  Nj  for  an  AR(1)  process  with  a  one-lag  temporal 
correlation  parameter  X=0.9999, 
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BIAS  OF  WHITE  NOISE  DRIVING 
VARIANCE  ESTIMATE 


X 


2 

Figure  17  Bias  of  for  the  time-averaged  parameter  estimators  of  an  AR(1) 
process  versus  the  one-lag  temporal  correlation  parameter  X  using 
N^IO. 
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A  parameter  estimators  of 
-empo-l  correlation  parameter  X 

using 
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0.1  10  -0.101 


-0.101 


0.5  10  -0.625 


-0.625 


0.6  10  -0.816 


0.7  10  -1.043 


-1.043 


0.8  10  -1.312 


-1.312 


0.85  10  -1.464 


-1.464 


0.9  10  -1.629 


-1.629 


0.99  10  -1.9603 


-1.9603 


0.9999  10  -1.9995 


-1.9995 


mean  Re  S(l) 
Burg  Algorith. 


-0.089682 


-0.099617 


-0.56279 


-0.61665 


-0.73757 


-0.816  -0.80760 


-0.94039 


-1.0317 


-1.1861 


-1.2994 


-1.3292  0.7225 


-1.4964 


-1.6133 


-1.8633 


-1.9470 


-1.9809 


-1.9974 


mean  Re  S(2) 
Burg  Algorith. 


0.012920 


0.0098845 


0.20679 


0.24397 


0.29377 


0.35263 


0.40322 


0.47984 


0.52919 


0.62764 


0.6045 


-1.4494  0.7225  0.70769 


0.6955 


0.7948 


0.89688 


0.96703 


3.9596 


3.9596 


2.8125 


2.8125 


2.2282 


2.2282 


1.5502 


1.5502 


0.8502 


0.8502 


0.53057 


0.53057 


0.26136 


0.26136 


3.1361e-3 


3.1361e-3 


3.1581 


3.8812 


2.249 


2.7597 


1.7904 


2.1835 


1.2639 


1.519 


0.70077 


0.83313 


0.44617 


0.52050 


0.22303 


0.25648 


4.4644x10-^ 


3.0977x10-- 


0.9997 


0.98197  4.7677e-7  1.9573x10-^ 


0.9997  0.99774  4.7677e-7  6.3050x10- 


Table  5a  Mean  of  estimated  parameters  for  a  Gaussian  AR(2)  process  using  the 
Burg  algorithm  with  X  and  Nj  as  parameters  and  Nj^=  10,000. 
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X 

eiTor  varS(l) 
Burg  Algorith. 

error  var  S(2) 
Burg  Algorith. 

error  var  6u 
Burg  Algorith 

0.1 

10 

1.1012x10-1 

1.1393x10-1 

1.9373 

100 

1.0019x10-2 

1.0303x10-2 

1.6370x10-1 

0.5 

10 

1.1239x10-2 

1.1166x10-1 

9.398x10-1 

100 

9.6568x10-3 

9.5378x10-3 

8.0525x10-2 

0.6 

10 

1.1145x10-1 

1.0946x10-1 

5.9132x10-1 

100 

9.1010x10-3 

9.0057x10-3 

5.0494x10-2 

0.7 

10 

1.0683x10-1 

1.0539x10-1 

2.8203x10-1 

100 

7.9729x10-3 

8.040x10-3 

2.4269x10-2 

0.8 

10 

9.9714x10-2 

9.6516x10-2 

8.5696x10-2 

100 

6.4081x10-3 

6.2261x10-3 

7.3920x10-3 

0.85 

10 

9.39520x10-2 

8.9116x10-2 

3.3739x10-2 

100 

5.2861x10-3 

5.2873x10*3 

2.8461x10-3 

0.9 

10 

7.992x10-2 

7.3954x10-2 

8.181x10-3 

100 

4.0572x10-3 

4.006x10-3 

7.026x10-4 

0.99 

10 

3.2103x10-2 

3.0368x10-2 

1.2169x10-5 

100 

8.405x10-4 

8.3077x10-4 

9.938x10-8 

0.9999 

10 

5.6797x10-3 

4.6069x10-3 

1.7209x10-111 

100 

5.0524xx10-3 

5.0442x10-5 

1.1187x10-13 

Table  5b  Error  variances  of  the  estimated  parameters  for  a  Gaussian  AR(2) 
process  using  the  Burg  algorithm  with  X  and  Nj  as  parameters  and 

Nr=  10,000. 
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X 

Nj 

sample  varS(l) 
Burg  Algorith. 

sample  var  S(2) 
Burg  Algorith. 

sample  var 

Burg  Algorith 

0.1 

10 

1.1000x10-^ 

1.1394x10*' 

1.2950 

100 

1.0018x10*2 

1.0304x10*2 

1.5756x10*' 

0.5 

10 

1.0851x10*' 

1.0978x10*' 

6.2231x10*' 

100 

9.5875x10*^ 

9.5023x10*2 

7.7746x10*2 

0.6 

10 

1.0530x10*' 

1.0508x10*' 

3.9965x10*' 

100 

9.0303x10*2 

8.9521x10*2 

4.8494x10*2 

0.7 

10 

9.6310x10*2 

9.7858x10*2 

2.0006x10*' 

100 

7.8448x10*2 

7.9375x10*2 

2.3301x10*2 

0.8 

10 

8.3871x10*2 

8.4242x10*2 

6.3381x10*2 

100 

6.2500x10*2 

6.0740x10*2 

7.1021x10*2 

0.85 

10 

7.5749x10*2 

7.5200x10*2 

2.6618x10*2 

100 

5.0711x10*2 

5.0684x10*2 

2.7448x10*2 

0.9 

10 

6.2334x10*2 

6.0849x10*2 

6.7118x10*2 

100 

3.8104x10*2 

3.7745x10*2 

6.7881x10*^ 

0.99 

10 

2.2685x10*2 

2.3444x10*2 

1.0406x10*5 

100 

6.6275x10*^ 

6.5986x10*^ 

9.7917x10*8 

0.9999 

10 

4.3379x10*2 

4.2963x10*2 

1.7008x10*"^ 

100 

4.6262x10*5 

4.6589x10*5 

8.8243x10*'^ 

Table  5c  Sample  variances  of  the  estimated  parameters  for  a  Gaussian  AR(2) 
process  using  the  Burg  algorithm  with  X  and  N-j-  as  parameters  and 

Nr= 10,000. 
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estimator.  The  estimates  were  obtained  using  batch  sizes  of  Nj=10  and  100. 
Again,  the  number  of  repeated  realizations  was  Njj=10,000.  Figs.  19,  20  and  21 

show  the  error  variances  of  ^[(1),  ^(2)  and  6^,  respectively.  Again,  a  significant 
decrease  in  the  error  variance  is  noted  at  high  values  of  temporal  correlation. 
For  ^u,  the  reduction  in  error  variance  using  Nj=100  is  approximately  12  orders 

of  magnitude  between  X,=0.1  and  0.9999. 

The  important  point  to  be  made  regarding  the  above  observations  is  that 

the  error  variance,  the  estimator  variance  and  the  bias  of  the  estimators  are  not 
only  dependent  upon  Nj,  but  also  on  process  correlation.  Furthermore,  the 

superiority  of  a  given  estimator  may  change  significantly  depending  upon  process 
correlation  especially  for  low  valuses  of  Nj. 
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Figure  19  Error  variance  of  ^i(l)  for  an  AR(2)  process  using  the  Burg  algorithm 


with  Nj=10  and  100. 


4  0.6  0.8  1.0 


Figure  20  Error  variance  of  ^(2)  for  an  AR(2)  process  using  the  Burg  algorithm 


with  Nj=10  and  100. 
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ERROR  VARIANCE  OF 
WHITE  NOISE  DRIVING  VARIANCE 


Figure  21  Error  variance  of  6$  for  an  AR(2)  process  using  the  Burg  algorithm 
with  Nj=10  and  100. 
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5.3.2  Performance  with  Non-Gaussian  SIRP  Processes 

In  this  subsection,  computed  simulation  results  of  the  error  variances  are 
presented  using  the  Burg  estimator  for  autoregressive  (AR)  time  series  models 
where  the  white  noise  driving  process  is  non-Gaussian.  Again,  the  order  of  the 
process  is  assumed  known.  The  non-Gaussian  processes  are  modeled  from  the 
general  class  known  as  spherically  invariant  random  processes  (SIRP).  In  [7],  the 
synthesis  of  such  processes  is  developed.  The  special  case  of  a  K-distributed  SIRP 
process  is  considered  here.  Detection  analyses  involving  K-distributed  clutter 
processes  are  presented  in  [8].  For  processes  consisting  of  in-phase  (real)  and 
quadrature  (imaginary)  components,  the  K-distributed  envelope  PDF  is  expressed 
as 


W  = 


2Va  Vai 
r(a)L  2  J 


a 


K„.,(>far) 


(0  <  r  <  oo) 


where  r(a)  is  the  Eulerian  Gamma  function,  K„(«)  is  the  modified  Bessel 

function  of  the  second  kind  with  order  a.  Here,  a  is  refered  to  as  the  shape 
parameter. 

AR(4)  processes  with  various  degrees  of  temporal  correlation  are 
synthesized  using  various  shape  parameters,  a.  The  sample  mean,  sample 
variance  and  error  variances  of  the  estimated  AR  parameters  and  white  noise 
driving  variances  are  listed  in  Tables  6a  through  6d  for  several  values  of  a  and 
X.  In  all  cases,  Nj=1000  time  samples  were  used  to  obtain  each  estimate  and 
Nr=1000  realizations  were  used  to  obtain  the  statistics. 

Examination  of  the  tables  reveals  that  using  N-j-=1000  time  samples,  good 
estimates  with  low  error  variances  are  obtained  for  the  a(k)  coefficients,  k=l  to 
4,  for  nearly  all  values  of  a  and  X.  Furthermore,  the  error  variances  for  the 
estimates  of  a(k)  are  insensitive  to  a.  The  same  point  cannot  be  made  for  the 
estimate  of  the  white  noise  driving  variance,  dy,  however.  Although,  the 
estimates  of  the  mean  values  of  are  quite  good  (using  Nj=l(X)0),  the  associated 
error  variances  are  quite  revealing.  In  Figs.  18  and  19,  we  plot  the  error 
variances  of  as  a  function  of  a  with  X,  as  a  parameter.  Each  curve  reveals  a 
significant  decrease  in  the  error  variance  with  increasing  a  (approximately  two 
orders  of  magnitude  between  a=0.1  and  10).  TTiis  result  is  explained  by  first 
noting  that  in  the  limit  as  a  ->  oo,  K-distributed  processes  approach  that  of  the 
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a 

sample 

mean 

sample 

variance 

error 

variance 

bias 

a(l) 

-0.32965 

-0.32885 

9.6365x10-4 

9.6442x10-4 

+8.000x10-4 

a(2) 

0.09963 

0.099624 

1.0361x10-3 

1.0361x10-3 

-6.000x10-6 

0.1 

a(3) 

-0.02967 

-0.030629 

1.0648x10-3 

1.0665x10-3 

-9.590x10-4 

a(4) 

0.0081 

0.009436 

1.0500x10-3 

1.0518x10-3 

+  1.336x10-3 

2 

3.6076 

3.7466 

169.19 

169.21 

0.139 

a(l) 

-0.32965 

-0.32946 

1.0004x10-3 

1.0004x10-3 

+1.900x10-4 

a(2) 

0.09963 

0.099045 

1.0895x10-3 

1.0899x10-3 

-5.850x10-4 

0.5 

a(3) 

-0.02967 

-0.02976 

1.0977x10-3 

1.0978x10-3 

-9.000x10-5 

a(4) 

0.0081 

0.008522 

9.9183x10-4 

9.9223x10-4 

+4.220x10-4 

2 

3.6076 

3.6551 

26.818 

26.820 

+4.750x10-2 

a(l) 

-0.32965 

-0.3302 

9.4419x10-4 

9.4479x10-4 

-5.500x10-4 

a(2) 

0.09963 

0.09920 

1.0639x10-3 

1.0644x10-3 

-4.300x10-4 

1.0 

a(3) 

-0.02967 

-0.028524 

1.1300x10-3 

1.1320x10-3 

+1.146x10-3 

a(4) 

0.0081 

0.007151 

9.7787x10-4 

9.7973x10-4 

-9.490x10-4 

2 

3.6076 

3.6619 

13.606 

13.609 

+5.430x10-2 

a(l) 

-0.32965 

-0.33011 

9.5535x10-4 

9.5572x10-4 

-4.600x10-4 

a(2) 

0.09963 

0.10060 

1.0742x10-3 

1.0752x10-3 

+9.700x10-4 

10 

a(3) 

-0.02967 

-0.029265 

1.0702x10-3 

1.0752x10-3 

+4.050x10-4 

BEII 

0.0081 

0.0078392 

1.0327x10-3 

1.0328x10-3 

-2.610x10-4 

2 

3.6076 

3.6188 

1.237 

1.2378 

+1.120x10-2 

a(l) 

-0.32965 

-0.32947 

1.0314x10-3 

1.0310x10-3 

+1.800x10-4 

a(2) 

0.09963 

0.098557 

1.0885x10-3 

1.0886x10-3 

-1.073x10-3 

oo 

a(3) 

-0.02967 

-0.028013 

1.1002x10-3 

1.1034x10-3 

+1.657x10-3 

a(4) 

0.0081 

0.0066742 

1.0106x10-3 

1.0124x10-3 

-1.426x10-3 

2 

3.6076 

3.5965 

1.3756x10-2 

1.3866x10-2 

-1.110x10-2 

Table  6a.  Performance  of  the  Burg  estimator  for  an  AR(4)  process  with  a  as  a 
parameter  and  A.=0.3,  Ny=1000  and  Nj^=1000. 
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a(l) 
a(2) 
1.0  I  a(3) 


-1.2934 


1.0513 


-0.63378 


0.2401 


1.2891 


-1.2934 


1.0513 


-0.63378 


0.2401 


1.2890 


-1.2934 


1.0513 


-0.63378 


0.2401 


1.2890 


-1.2934 


1.0513 


-0.63378 


0.2401 


1.2890 


-1.2934 


1.0513 


-0.63378 


0.2401 


1.2890 


sample 

mean 

sample 

variance 

error 

variance 

bias 

-1.2914 

9.4170x10-4 

9.456x10-4 

+2.000x10-3 

1.0486 

2.0586x10-3 

2.066x10-3 

-2.700x10-3 

-0.6312 

2.0904x10-3 

2.099x10-3 

+2.580x10-3 

0.2388 

9.0343x10-4 

-1.300x10-3 

1.2810 

3.137x10-2 

3.144x10-2 

-8.100x10-3 

-1.2911 

9.3737x10-4 

9.4255x10-4 

+2.300x10-3 

1.0492 

2.1866x10-3 

2.1910x10-3 

-2.100x10-3 

-0.63199 

2.2355x10-3 

2.2390x10-3 

+1.790x10-3 

0.23968 

9.0611x10-4 

9.0646x10-4 

-4.200x10-4 

1.2807 

8.1152x10-3 

8.1839x10-3 

-8.300x10-3 

-1.2925 

9.0805x10-4 

9.0888x10-4 

+9.000x10-4 

1.0485 

2.0950x10-3 

2.1027x10-3 

-2.800x10-3 

-0.6299 

2.1903x10-3 

2.2059x10-3 

+3.880x10-3 

0.2377 

9.3840x10-4 

9.4483x10-4 

-2.400x10-3 

1.2841 

4.9590x10-3 

4.9826x10-3 

-4.900x10-3 

-1.2917 

9.3806x10-4 

9.4155x10-4 

+1.700x10-3 

1.0484 

2.2654x10-3 

2.2762x10-3 

-2.900x10-3 

-0.63092 

2.3454x10-3 

2.3551x10-3 

+2.860x10-3 

0.2383 

9.8906x10-4 

9.9259x10-4 

-1.800x10-3 

1.2830 

2.0194x10-3 

2.0554x10-3 

-6.000x10-3 

-1.2923 

9.0910x10-4 

9.0947x10-4 

+1.100x10-3 

1.0497 

2.0942x10-3 

2.0951x10-3 

-1.600x10-3 

-0.63156 

2.1302x10-3 

2.1340x10-3 

+2.220x10-3 

0.23923 

9.2802x10-4 

9.2785x10-4 

-8.700x10-4 

1.2834 

1.6744x10-3 

1.7038x10-3 

-5.600x10-3 

Table  6b.  Performance  of  the  Burg  estimator  for  an  AR(4)  process  with  a  as  a 
parameter  and  X=0J,  Nj=1000  and  Nj^=1000. 
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a 

sample 

mean 

sample 

variance 

error 

variance 

bias 

a(l) 

-2.6978 

-2.6949 

5.9992x10-4 

6.0853x10-4 

+2.900x10-3 

a(2) 

3.3081 

3.3013 

3.2060x10^ 

3.2524x10-3 

-6.800x10-3 

0.1 

a(3) 

-2.1852 

-2.1788 

3.1615x10-3 

3.2027x10-3 

+6.400x10-3 

a(4) 

0.6561 

0.65389 

5.8826x10-4 

5.9325x10-4 

-2.210x10-3 

2 

0.069747 

0.067509 

5.8782x10-2 

5.8787x10-2 

-2.238x10-3 

a(l) 

-2.6978 

-2.6939 

5.6190x10-4 

5.7698x10-4 

+3.900x10-3 

a(2) 

3.3081 

3.2985 

3.1312x10-3 

3.2225x10-3 

-9.600x10-3 

0.5 

a(3) 

-2.1852 

-2.1753 

3.1650x10-3 

3.2641x10-3 

+9.90x10-3 

a(4) 

0.6561 

0.6519 

5.9119x10-4 

6.0818x10-4 

-4.200x10-3 

2 

Ou 

0.069747 

0.070988 

9.9154x10-3 

9.9169x10-3 

+1.241x10-3 

a(l) 

-2.6978 

-2.6945 

5.5612x10-4 

5.6765x10-4 

+3.300x10-3 

a(2) 

3.3081 

3.3006 

3.0109x10-3 

3.0704x10-3 

-7.500x10-3 

1.0 

a(3) 

-2.1852 

-2.1778 

3.0589x10-3 

3.1154x10-3 

+7.400x10-3 

a(4) 

0.6561 

0.65298 

5.6884x10-4 

5.7859x10-4 

-3.120x10-3 

2 

0.069747 

0.07194 

5.1243x10-3 

5.1291x103 

+2.193x10-3 

a(l) 

-2.6978 

-2.6953 

5.7912x10-4 

5.8551x10-4 

+2.500x10-3 

a(2) 

3.3081 

3.3014 

3.0688x10-3 

3.1151x10-3 

-6.700x10-3 

10 

a(3) 

-2.1852 

-2.1781 

2.9719x103 

3.0236x10-3 

+7.100x10-3 

a(4) 

0.6561 

0.65319 

5.3793x10-4 

5.4664x10-4 

-2.910x10-3 

2 

0.069747 

0.069316 

4.9057x10-4 

4.9076x10-4 

-4.310x10-4 

a(l) 

-2.6978 

-2.695 

5.5435x10-4 

5.6151x10-4 

+2.800x10-3 

3.3081 

3.3008 

2.9622x10-3 

3.0152x10-3 

-7.300x10-3 

oo 

a(3) 

-2.1852 

-2.1773 

2.9584x10-3 

3.0199x10-3 

+7.900x10-3 

a(4) 

0.6561 

0.6526 

5.6196x10-4 

5.7378x10-4 

-3.500x10-3 

2 

0.069747 

0.069499 

4.9018x10-6 

4.9591x10-6 

-2.480x10-4 

Table  6c.  Performance  of  the  Burg  estimator  for  an  AR(4)  process  with  a  as  a 
parameter  and  A.=0.9,  Nj=1000  and  Nj^=1000. 
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a 

sample 

mean 

sample 

variance 

error 

variance 

bias 

a(l) 

-3.8729 

-3.8696 

2.7532x10-6 

3.8603x10-6 

+3.30x10-3 

a(2) 

5.7379 

5.7280 

2.3402x10-4 

3.3179x10-4 

-9.900x10-3 

0.1 

a(3) 

-3.8542 

-3.8442 

2.3299x10-4 

3.3238x10-4 

+1.100x10-2 

a(4) 

0.99033 

0.98695 

2.7189x10-6 

3.8638x10-6 

-3.380x10-3 

2 

3.4907x10-6 

2.8519x10-6 

7.5142x10-" 

7.555x10-11 

-6.388x10-2 

a(l) 

-3.8729 

-3.8698 

2.6899x10-6 

3.6670x10-6 

+3.100x10-3 

a(2) 

5.7379 

5.7286 

2.2970x10-4 

3.1683x10-4 

-9.300x10-3 

0.5 

a(3) 

-3.8542 

-3.8448 

2.2980x10-4 

3.1921x10-4 

+9.400x10-3 

a(4) 

0.99033 

0.98711 

2.6922x10-6 

3.7300x10-6 

-3.220x10-3 

2 

Ou 

3.4907x10-6 

3.5639x10-6 

2.6984x10-" 

2.6990x10-11 

+7.320x10-^ 

a(l) 

-3.8729 

-3.8698 

2.6217x10-6 

3.5597x10-6 

+3.100x10-3 

a(2) 

5.7379 

5.7287 

2.2341x10-4 

3.0717x10-4 

-9.200x10-3 

1.0 

a(3) 

-3.8542 

-3.8449 

2.2297x10-4 

3.0913x10-4 

+9.300x10-3 

a(4) 

0.99033 

0.98717 

2.6076x10-6 

3.6103x10-6 

-3.160x10-3 

2 

3.4907x10-6 

3.1827x10-6 

9.7524x10-12 

9.8023x10-12 

-3.080x10-2 

a(l) 

-3.8729 

-3.8697 

2.6989x10-6 

3.7185x10-6 

+3.200x10-3 

a(2) 

5.7379 

5.7284 

2.2932x10-4 

3.2023x10-4 

-9.500x10-3 

10 

a(3) 

-3.8542 

-3.8445 

2.2859x10-4 

3.2206x10-4 

+9.700x10-3 

a(4) 

0.99033 

0.98703 

2.6732x10-6 

3.7627x10-6 

-3.300x10-3 

2 

3.4907x10-6 

3.5003x10-6 

1.2889x10-12 

1.2890x10-12 

+9.600x10-9 

a(l) 

-3.8729 

-3.8698 

2.5457x10-6 

3.4981x10-6 

+3.100x10-3 

a(2) 

5.7379 

5.7288 

2.1643x104 

2.9991x10-4 

-9.100x10-3 

oo 

a(3) 

-3.8542 

-3.8449 

2.1575x10-4 

3.0044x10-4 

+9.300x10-3 

a(4) 

0.99033 

0.98717 

2.5218x10-6 

3.5166x10-6 

-3.160x10-3 

2 

3.4907x10-6 

3.4796x10-6 

1.1347x10-14 

1.1458x10-14 

-1.110x10-8 

Table  6d.  Performance  of  the  Burg  estimator  for  an  AR(4)  process  with  a  as  a 
parameter  and  A,=0.99,  Nj=1000  and  N(^=1000. 
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Gaussian.  Thus,  for  small  values  of  a,  we  have  processes  which  depart  most 
significantly  from  Gaussianity.  For  these  processes,  the  tails  of  the  PDF 
distribution  are  highest.  Thus,  for  a  fixed  sample  size  N^,  the  uncertainty  in  the 

estimate  of  the  variance  is  expected  to  increase  as  a  decreases.  In  Fig.  19,  we 
show  an  expanded  view  of  the  upper  three  curves  from  Fig.  18.  In  this  figure,  we 
note  a  two  order  of  magnitude  decrease  in  the  error  variance  over  the  range  of  a 
values  from  a=0.1  to  10.  From  the  tables,  we  observe  that  these  curves  continue 
to  decrease  by  nearly  two  additional  orders  of  magnitude  as  a  — >  «».  We  also 
note  the  drastic  reduction  in  the  error  variance  with  increasing  temporal 
correlation  expressed  by  the  one-lag  temporal  correlation  parameter  X. 


a 

2 

Fig. 22  The  error  variance  of  the  estimate  versus  the  shape  parameter  a  with 
the  one-lag  temporal  correlation  parameter  X  fixed. 
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5.4  Computed  Error  Variance  of  Multichannel  AR  Model  Parameters 


In  this  section,  we  present  the  computed  results  for  the  means  and  error 

variances  associated  with  the  multichannel  estimates  of  the  JxJ  matrix  coefficients 
A(k)  and  [Zf]^  defined  by  eqs.(2.1)  and  (2.6),  respectively.  The  Strand-Nuttall 

algorithm  is  used  in  all  the  cases  presented  here.  Specific  consideration  is  given 
to  the  performance  of  the  estimator  not  only  as  a  function  of  the  temporal 
correlation  on  each  channel,  but  also  the  cross-channel  correlation.  We  consider 

multichannel  AR(2)  processes  with  J=2  channels  and  various  values  of  the  one-lag 
temporal  correlation  parameter  X,jj  on  channel  j  and  the  cross-channel  correlation 

parameter,  Ip-I.  In  this  case,  we  have  parameters 


and 


A(k)  = 


aii(k)ai2(k) 
a2i(k)  a22(k) 


^11  ^12 
.^21  ^22. 


k=l,2 


(5.4) 

(5.5) 


Table  7  contains  the  computed  results  for  the  estimates  Sii(l),  ^11 

and  ^12-  The  number  of  time  samples  used  to  obtain  each  estimate  was  Ny=100 
while  the  statistics  were  computed  with  Nr=10(X)  realizations.  In  Figs.  22  and 

23,  we  plot  the  error  variances  of  the  matrix  element  estimates  ^^(l)  and  tjj.  , 

Similar  results  were  obtained  for  the  estimates  of  the  other  elements. 

Consistent  with  the  results  noted  for  single  channel  processes,  the  error 
variances  associated  with  the  estimates  of  these  coefficients  decrease  with 
increasing  temporal  correlation.  However,  we  now  observe  in  Figs  22  and  23 
that  the  error  variance  associated  with  the  estimates  Sii(l)  and  ^12(1)  increases 

with  increasing  cross-channel  correlation.  In  Figs. 24  and  25,  however,  we 
observe  that  the  error  variance  associated  with  the  estimates  of  the  white  driving 
noise  covariance  matrix  elements  ilj,  and  t.^2  independent  of  the  cross¬ 
channel  correlation.  This  result  is  consistent  with  the  presentation  made  in  section 
3.2  that  error  variances  of  cross-correlation  function  estimates  are  independent  of 
the  cross-channel  correlation  for  the  special  case  of  wide-sense  jointly  stationary 
narrowband  bandpass  processes  as  considered  here. 
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^11 

^2 

^12 

'P12’ 

0.1 

0.1 

0.01 

0 

TRUE 

VALUE 

MEAN  OF 

ESTIMATE 

SAMPLE 

VARIANCE 

ERROR 

VARIANCE 

-0.101 

-0.099721 

1.0263x10-2 

0 

-2.1981x10-3 

1.0529x10-2 

1.0525x10-2 

3.9596 

3.8077 

1.4491x10-^ 

1.6783x10-* 

0 

HEBSSIBO 

1.4529x10-* 

1.4556x10-* 

-0.625 

-0.61265 

9.7337x10-3 

9.8814x10-3 

0 

4.3440x10-^ 

9.3721x10-3 

9.3647x10-3 

2.8125 

2.7058 

7.4709x10-2 

8.6021x10-2 

0 

1.4185x10-3 

7.7050x10-2 

7.7009x10-2 

-1.043 

-1.0209 

8.6619x10-3 

9.1400x10-3 

0 

8.3557x10-3 

8.3803x10-3 

1.5502 

1.4824 

2.3082x10-2 

2.7649x10-2 

0 

2.5259x10-3 

2.4681x10-2 

2.4676x10-2 

-1.6290 

-1.5929 

4.7884x10-3 

6.0843x10-3 

0 

-6.4223x10-4 

4.8277x10-3 

4.8239x10-3 

0.26136 

0.25166 

6.5474x10-4 

7.4833x10-4 

0 

1.9544x10-4 

6.5406x10-4 

6.5357x10-4 

-1.9603 

-1.9180 

1.8241x10-3 

3.6105x10-3 

0 

1.8338x10-5 

2.2002x10-3 

2.1980x10-3 

3.1361x10-^ 

3.0368x10-3 

9.3876x10-8 

1.0365x10-2 

0 

-6.4657x10-6 

9.8708x10-8 

9.8692x10-8 

Table  7a  Tabulated  values  of  the  mean  and  variances  for  the  multichannel 
coefficient  estimates  Sn(l),  ^12(1)  »  ^ii(0  and  ^12(1)  for  AR(2) 
processes  with  specified  temporal  and  cross-channel  correlation 
using  the  Strand-Nuttall  estimator. 
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A.JJ 

^2 

^12 

'P12’ 

0.1 

0.1 

0.1 

0.7 

0.5  0.5  0.5  0.7 


0.7  0.7  0.7  0.7 


0.9  0.9  0.9  0.7 


0.99  0.99  0.99  0.7 


TRUE 

MEAN  OF 

SAMPLE 

ERROR 

VALUE 

ESTIMATE 

VARIANCE 

VARIANCE 

-0.101 

-0.10209 

1.9722x10-2 

1.9706x10-2 

0 

5.6071x10-3 

2.0244x10-2 

2.0256x10-2 

3.9596 

3.8016 

1.4683x10-’ 

1.7166x10-’ 

2.7717 

2.6534 

1.4752x10-’ 

1.6142x10-’ 

-0.625 

-0.61097 

1.9408x10-2 

1.9587x10-2 

0 

2.0894x10-3 

1.9262x10-2 

1.9249x10-2 

2.8125 

2.7092 

7.8848x10-2 

8.9437x10-2 

1.9687 

1.9026 

7.3984x10-2 

7.8289x10-2 

-1.043 

-1.0235 

1.6809x10-2 

1.7196x10-2 

2.557x10-'^ 

1.5809x10-3 

1.7263x10-2 

2.6020x10-2 

1.5502 

1.4890 

2.3195x10-2 

2.6915x10-2 

1.0851 

1.0405 

2.4046x10-2 

2.6020x10-2 

-1.6290 

-1.5993 

8.7113x10-3 

9.9725x10-3 

-2.1458x10*'^ 

-7.7136x10-^ 

8.9255x10-3 

8.9177x10-3 

0.26136 

0.25214 

6.1706x10-^ 

7.0158x10-^ 

0.18295 

0.17707 

6.4987x10-^ 

6.8409x10-^ 

-1.9603 

-1.9168 

4.1238x10-3 

6.0096x10-3 

5.9009x10-'^ 

-7.0187x10-^ 

4.3465x10-3 

4.3438x10-3 

3.1361x10-3 

3.0342x10-3 

1.0060x10-2 

1.1090x10-2 

2.195x10-3 

2.1211x10-3 

9.8467x10-8 

1.0388x10-2 

Table  7b  Tabulated  values  of  the  mean  and  variances  for  the  multichannel 
coefficient  estimates  ^12(1) »  ^ii(0  and  AR(2)  processes  with 

specified  temporal  and  cross-channel  correlation  using  the  Strand-Nuttall 
estimator. 
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^11 

^2 

^12 

'P12' 

TRUE 

VALUE 

0.1 

0.1 

0.1 

0.99 

an 

-0.101 

^12 

4.8161x10'^ 

3.9596 

^12 

3.92 

0.5 

0.5 

0.5 

0.99 

ail 

-0.625 

ai2 

0 

2^11 

2.8125 

^12 

2.7844 

0.7 

0.7 

0.7 

0.99 

an 

-1.043 

ai2 

-1.503x10-5 

^11 

1.5502 

^12 

1.5347 

0.9 

0.9 

0.9 

0.99 

an 

-1.6290 

ai2 

1.0300x10-5 

2ii 

0.26136 

^12 

0.25873 

0.99 

0.99 

0.99 

0.99 

an 

-1.9603 

ai2 

-7.1755x10-5 

2^11 

3.1361x10-^ 

^12 

3.067x10-3 

Table  7c  Tabulated  values  of  the  mear 
coefficient  estimates  Sj/l),  ^12(0  ’ 
specified  temporal  and  cross-channel  c 
estimator. 
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MEAN  OF 

SAMPLE 

ERROR 

ESTIMATE 

VARIANCE 

VARIANCE 

-0.1027 

4.9512x10-^ 

4.9503x10-* 

3.1909x10-3 

4.9880x10-* 

4.9871x10-* 

3.8008 

1.5615x10-* 

1.8122x10* 

3.7651 

1.5696x10* 

1.8080x10* 

0.6230 

5.1454x10-* 

5.1421x10-* 

1.0188x10-2 

5.1596x10-* 

5.1579x10-* 

2.6910 

7.7955x10-2 

9.2635x10-2 

2.6650 

7.7097x10-2 

9.1265x10-2 

-1.0390 

4.2740x10-* 

4.2803x10  * 

1.6488x10-2 

4.2911x10-* 

4.3012x10-* 

1.4899 

2.3652x10-2 

2.7263x10-2 

1.4759 

2.3821x10-2 

2.7255x10-2 

-1.5976 

2.3052x10-* 

2.3128x10-* 

2.5663x10-3 

2.3220x10-* 

2.3197x10-* 

0.25186 

6.6874x10-^ 

7.5798x10-^ 

0.24928 

6.6927x10-^ 

7.5788x10-^ 

-1.9143 

1.1510x10* 

1.1712x10-* 

-3.7304x10-3 

1.1472x10-* 

1.1462x10-* 

2.9814x10-3 

1.0885x10-2 

1.2240x10-2 

2.9483x10-3 

1.1691x10-2 

1.3090x10-2 

and  variances  for  the  multichannel 
and  ^12(1)  for  AR(2)  processes  with 
)rrelation  using  the  Strand-Nuttall 


ERROR  VARIANCE 
FOR  all(l)  ESTIMATE 


X 


Figure  24  Error  variance  for  the  estimate  coefficient  versus  the  one-lag 
temporal  correlation  parameter  using  the  Strand-Nuttall 

algorithm  with  order  2,  Nj=100  time  samples,  and  Ipj2l  as  a 
parameter. 
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ERROR  VARIANCE 


Figure  25  Error  variance  for  the  estimate  coefficient  versus  the  one-lag 
temporal  correlation  parameter  A,=A,,j=A,,2  using  the  Strand-Nuttall 
algorithm  with  order  2,  Nj=100  time  samples,  and  lp,2l  as  a 

parameter. 
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6.0  SUMMARY 


In  this  study,  the  performance  of  time-averaged  estimators  is  considered 
not  only  as  a  function  of  the  observation  window  size  of  the  data,  but  also  in 
terms  of  the  pertinent  correlation  parameters  of  the  underlying  observation 
processes.  Specifically,  the  error  variances  and  sample  variances  of  both 
correlation  function  and  parameter  estimators  are  considered  for  Gaussian  as  well 
as  non-Gaussian  processes.  In  addition,  multichannel  processes  are  analyzed. 
Analytic  expressions  are  derived  for  the  variance  of  the  time-averaged  complex 
auto-  and  cross-correlation  functions.  The  unique  aspect  of  this  development  is 
the  determination  of  the  functional  dependence  of  these  expressions  in  terms  of 
process  correlation  parameters.  Specifically,  the  variance  of  the  time-averaged 
cross-correlation  function  is  shown  to  depend  upon  the  temporal  correlation  and 
variance  of  each  process,  as  well  as  (in  the  more  general  case  of  processes  with 
unconstrained  Gaussian  quadrature  components  [5])  the  cross-correlation 
coefficient.  These  expressions  provide  a  performance  measure  which  can  be  used 
to  specify  the  window  size  of  the  observation  interval  required  to  achieve  a 
specific  value  of  this  variance. 
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APPENDIX  A 


In  this  Appendix,  we  derive  the  result  expressed  in  eq.(29b)  of  the  text. 
From  eqs(l  1) 

E[(l)(n,l)]  =  Rii(l)  (A.la) 

and 

E[(l)*(n-k,l)|  =  R‘i(l)  (A.lb) 

SO  that  from  eq(2d) 

C,,,^,(k,l)  =  R^^(k,l)  -  lRii(l)P.  (A.2) 

And  SO,  eq(3c)  becomes 

Vb.O.Nt)  =  if  I  [l  -  [R„(k.l)  -IR„(I)P].  (A.3) 

T  k=-(N.j.-ll|-l)‘-  -• 

We  now  consider, 

R^^(k,l)  =  El(|)(n,l)(t.*(n  -  k,l)l  (A.4a) 

=  E[xj(n)xi  (n  -  l)xi  (n  -  k)xj(n  - 1  -  k)].  (A.4b) 

For  processes  with  zero-mean,  jointly  stationary  Gaussian  quadrature 
components  Xjj(n)  and  XiQ(n),  eq(A.4b)  can  be  expressed  as  [see  Appendix  B] 

R^$(k,l)  =  I  Rii(l)|2  +  I  Ri;(k)P  +  Fii(l,k)  (A.5) 

where 

Fii(l,k)  =  ElXi(n)Xi(n  - 1  -  k)JE[x*(n  -  l)x*(n  -  k))  (A.6a) 

=  {  r"(I  +  k)  -  R^’^d  +  k))  { r"(I  -  k)  -  R^'^d  -  k) } 

+  { R^'d  +  k)  +  rIPcI  +  k) )  { R^'d  -  k)  +  R'iPd  -  k) } 

-  j  [  R  "(I  +  k)  -  RP'^d  +  k) )  {  rP'(I  -  k)  +  R'iPd  -  k) ) 

+  j  (Rud-k)-Rp‘^d-k))lRa'd  +  k)  +  R|Pd  +  k))-  (A.6b) 

Using  eq(A.5)  in  (A.3),  we  have 
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(A,7) 


VehC-Nt)  =  if  I  1  -  [lRil(k)P  +  Fi,(l,k)]. 

“  "Tk=-(N^-ll|-l)'-  tJ 

By  examination  of  eq(A.6b),  we  note  that  the  imaginary  terms  in  eq(A.7) 
sum  to  zero.  This  can  be  seen  by  first  noting  that  Fjj(l,k)  is  real  for  k=0.  We  also 

note  that  imaginary  terms  evaluated  with  negative  values  of  k  serve  to  cancel  the 

corresponding  imaginary  terms  for  positive  values  of  k.  And  so,  only  the  real 
part  of  the  function  Fjj(l,k)  contributes  to  the  Vg..(l,Nj)  function.  Therefore, 

N.j,-ll|-l 

VBiiO.NT)  =  if  I  I  [lRii(k)P  +  RelFii(l,k))]  (A.8) 

Tk=-(N^-ll|-l)L  fJ 

which  is  eq(4)  in  the  text. 
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APPENDIX  B 

In  this  appendix,  we  derive  eq(A.5)  of  Appendix  A.  Consider  eq(A.4b) 
expressed  as 

R(])(j,(k,l)  =  E[xi(n)xi  (n  -  l)xi  (n  -  k)xi(n  - 1  -  k)].  (B.  1 ) 

In  the  special  case  where  the  process  xj(n)  is  Gaussian,  then  [6] 

R^^(k,l)  =  E[Xi(n)x- (n  -  l)]E[x*(n  -  k)xi(n  - 1  -  k)] 

+  E[Xj(n)xi  (n  -  k)]E[xi  (n  -  l)Xj(n  - 1  -  k)].  (B.2) 

However,  we  do  not  wish  to  constrain  this  discussion  to  this  restrictive  case. 
Rather,  we  wish  to  consider  the  more  general  case  of  a  process  Xj(n)  with  jointly 

Gaussian  quadrature  components.  We  therefore  consider 

Xi(n)  =  Xji(n)  +  j  XiQ(n)  (B.3) 

where  the  processes  Xij(n)  and  XjQ(n)  are  jointly  Gaussian.  Using  eq(B.3)  in 
(B.l),  we  obtain 


R<j,<f,(k)  =  E  { [Xii(n)  +  j  XiQ(n)][Xii(n-l)  -  j  XjQCn-l)] 

•  [Xji(n-k)  -  j  XjQ(n-k)][Xji(n-l-k)  +  j  XiQ(n-l-k)] }  (B.4a) 

=  E  { [Xii(n)Xii(n-l)+XjQ(n)xiQ(n-l)+jXjQ(n)Xii(n-l)'jxii(n)xjQ(n-l)] 

.  [xii(n-k)Xij(n-l-k)+XjQ(n-k)xjQ(n-l-k)+jXji(n-k)XjQ(n-l-k) 

-jXiQ(n-k)Xji(n-l-k)] }  (B.4b) 

=  E[Xii(n)Xii(n-l)Xii(n-k)Xji(n-l-k)]  +  E[XjQ(n)XiQ(n-l)Xji(n-k)Xii(n-l-k)] 

+  E[Xji(n)xji(n-l)xiQ(n-k)xjQ(n-l-k)]  +  E[XjQ(n)xjQ(n-l)xiQ(n-k)XjQ(n-l-k)] 

-  E[XiQ(n)xji(n-i)Xii(n-k)XiQ(n-l-k)]  +  E[Xi,(n)XiQ(n-l)Xii(n-k)XiQ(n-l-k)] 

+  E[XiQ(n)Xii(n-l)XiQ(n-k)Xii(n-l-k)]  -  E[Xi,(n)XiQ(n-l)XiQ(n-k)Xji(n-l-k)] 
+jE[Xji(n)Xii(n-l)Xji(n-k)XjQ(n-l-k)]  -  jE[Xii(n)Xji(n-l)XiQ(n-k)Xji(n-l-k)] 
+jE[XiQ(n)xjQ(n-l)Xij(n-k)xjQ(n-l-k)]  -  jE[XiQ(n)XiQ(n-l)xjQ(n-k)Xii(n-l-k)] 
+jE[XjQ(n)Xji(n-l)Xji(n-k)Xii(n-l-k)]  +  jE[xjQ(n)Xii(n-l)xiQ(n-k)XiQ(n-l-k)] 
-jE[xji(n)XjQ(n-l)Xji(n-k)Xii(n-l-k)]  -  j  E[Xji(n)xiQ(n-l)XjQ(n-k)XjQ(n-l-k)] 

(B.4c) 

For  Gaussian,  zero-mean  quadrature  components,  eq(B.4c)  can  be 
expressed  as 
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RH(k)  =  [Ra  +  [Ru  (k)l^  +  Rad  +  k)Ra  (k  - 1) 

+  R§^(l)Ra  d)  +  [Ra '(k)l^  +  Ra  V  +  k)Ra '(k  - 1) 

+  Ra  d)R^d)  +  (RS^(k)]2  +  R'^d  +  k)R!P(k  - 1) 

+  IRa^d)]^  +  [R?^(k)l^  +  Ra^d  +  k)R^(k  - 1) 

-  R^^DRiPd)  -  R?(k)R!P(k)  -  Ra^d  +  k)Ra (k  - 1) 

+  [RlPd)]2  +  Ru(k)R?^(k)  +  R^Pd  +  k)RP'(k  - 1) 

+  [Ra '(1)1^  +  Ra^(k)Ra  (k)  +  Ra 'd  +  k)R!P(k  - 1) 

-  R!Pd)Ra 'd)  -  Rf  (k)Ra '(k)  -  Ra  d  +  k)Ra^(k  - 1) 

+j  { Ra  d)R!?d)  +  r" (k)R!P(k)  +  rIPcI  +  k)Ra  (k  - 1) ) 

-j  (Rad)Ra‘d)  +  RP(k)Ra(k)  +  Rad  +k)R'^(k  - 1)} 

+j  { Rp^CDR^^d)  +  Ra '(k)Ra^(k)  +  Ra^d  +  k)Ra '(k  - 1) ) 

-j  {RP'^d)RP'd)  +  RP^(k)Rp'(k)  +  Ra'd  +  k)Ra^(k  - 1)} 

+j  { Ra 'd)Ra  d)  +  Ra '(k)Ra  (k)  +  Rp 'd  +  k)RP(k  - 1) } 

+j  {Ra'd)Ra^d)  +  Ra^RiPdc)  +  Rp^d  +  k)R|P(k  - 1)} 

-j  { RlPd)RP(l)  +  Ra’(k)Rp‘(k)  +  Ra  d  +  k)Rp'(k  - 1) ) 

-j  { R  iPd)Rp^d)  +  r’^ (k)Rp^(k)  +  Ra^d  +  k)Rp^(k  - 1) ) 

(B.5) 

where  we  note  that  the  first  two  terms  in  each  parenthesis  for  the  imaginary 
terms  cancel.  Since 

Rd)  =  [Ra  d)  +  Ra^d)]  +  j  [Rp'd)  -  Ra^d)]  (B-6) 

then 

lRd)l^  =  (Ra  d)l^  +  2Ra  (DRp^d)  +  (Rp^d)]^ 

+  (Rp'(k)}^  -  2RP'd)RL‘^d)  +  [R^ (1)1^  (B.7) 

and  similarly  for  lR(k)P  so  that 

R(l,(„(k)  =  iRiid)!^  +  lRii(k)|2  +  Fii(l,k)  (B.8) 
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where 

Fii(l,k)  =  {Rad  +  k)  -  R^'^d  +  k))(  Rad  -  k)  -  R^'^d  -  k)} 

+  { R§'d  +  k)  +  RiPd  +  k) }  {  R^'d  -  k)  +  RlPd  -  k) ) 

-  j  {  r" (I  +  k)  -  R^'^d  +  k) )  (  R?'d  -  k)  +  RiPd  -  k) } 

+  j  {  Ra  (I  -  k)  -  R^d  -  k) )  {  R^'d  +  k)  +  RiPd  +  k) } 
as  noted  in  eq(A.6b)  of  Appendix  A. 
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Rome  Laboratory  plans  and  executes  an  interdisciplinary  program  in  re¬ 
search,  development,  test,  and  technology  transition  in  support  of  Air 
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Force  Command,  Control,  Communicatims  and  Intelligence  (Cl)  activities 
for  all  Air  Force  platforms.  It  also  executes  selected  acquisition  programs 
in  several  areas  of  expertise.  Technical  and  engineering  support  within 
areas  of  competence  is  provided  to  ESD  Program  Offices  (POs)  and  other 
ESD  elements  to  perform  effective  acquisition  of  C  I  systems.  M  addition, 
Rome  Laboratory's  technology  supports  other  AFSC  Product  Divisions,  the 
Air  Force  user  community,  and  other  DOD  and  non-DOD  agencies.  Rome 
Laboratory  maintains  technical  competence  and  research  programs  in  areas 
including,  but  not  limited  to,  communications,  command  and  control,  battle 
management,  intelligence  information  processing,  computational  sciences 
and  software  producibility,  wide  area  surveillance/sensors,  signal  proces¬ 
sing,  solid  state  sciences,  photonics,  electromagnetic  technology,  super¬ 
conductivity,  and  electronic  reliability/maintainability  and  testability. 


